module beams use const_and_precisions, only : wp_, one implicit none contains subroutine read_beam0(file_beam,fghz,x00,y00,z00, & wcsi,weta,rcicsi,rcieta,phiw,phir,unit) use const_and_precisions, only : pi,vc=>ccgs_ use utils, only : get_free_unit implicit none ! arguments character(len=*), intent(in) :: file_beam real(wp_), intent(out) :: fGHz,x00,y00,z00 real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw integer, intent(in), optional :: unit ! local variables integer :: u real(wp_) :: ak0,zrcsi,zreta,w0csi,w0eta,d0csi,d0eta if (present(unit)) then u=unit else u = get_free_unit() end if open(unit=u,file=trim(file_beam),status='OLD',action='READ') ! fghz wave frequency (GHz) read(u,*) fGHz ! x00,y00,z00 coordinates of launching point in cm read(u,*) x00, y00, z00 ! beams parameters in local reference system ! w0 -> waist, d0 -> waist distance from launching point ! phiw angle of spot ellipse read(u,*) w0csi,w0eta,d0csi,d0eta,phiw close(u) ak0=2.0e9_wp_*pi*fghz/vc zrcsi=0.5_wp_*ak0*w0csi**2 zreta=0.5_wp_*ak0*w0eta**2 wcsi=w0csi*sqrt(1.0_wp_+(d0csi/zrcsi)**2) weta=w0eta*sqrt(1.0_wp_+(d0eta/zreta)**2) rcicsi=-d0csi/(d0csi**2+zrcsi**2) rcieta=-d0eta/(d0eta**2+zreta**2) phir=phiw end subroutine read_beam0 subroutine read_beam1(file_beam,alpha0,beta0,fghz,x00,y00,z00, & wcsi,weta,rcicsi,rcieta,phiw,phir,unit) use const_and_precisions, only : pi,vc=>ccgs_ use simplespline, only : spli, difcs use utils, only : get_free_unit,locate implicit none ! arguments character(len=*), intent(in) :: file_beam real(wp_), intent(in) :: alpha0 real(wp_), intent(out) :: fghz,x00,y00,z00,beta0 real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw integer, intent(in), optional :: unit ! local variables integer :: u,iopt,ier,nisteer,i,k,ii real(wp_) :: steer,dal real(wp_), dimension(:), allocatable :: alphastv,betastv,x00v,y00v, & z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v, & cbeta,cx0,cy0,cz0,cwaist1,cwaist2, & crci1,crci2,cphi1,cphi2 if (present(unit)) then u=unit else u = get_free_unit() end if open(unit=u,file=file_beam,status='OLD',action='READ') read(u,*) fghz read(u,*) nisteer allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), & waist2v(nisteer),rci1v(nisteer),rci2v(nisteer), & phi1v(nisteer),phi2v(nisteer),x00v(nisteer), & y00v(nisteer),z00v(nisteer),cbeta(4*nisteer), & cx0(4*nisteer),cy0(4*nisteer),cz0(4*nisteer), & cwaist1(4*nisteer),cwaist2(4*nisteer),crci1(4*nisteer), & crci2(4*nisteer),cphi1(4*nisteer),cphi2(4*nisteer)) do i=1,nisteer read(u,*) steer,alphastv(i),betastv(i),x00v(i),y00v(i),z00v(i), & waist1v(i),waist2v(i),rci1v(i),rci2v(i),phi1v(i),phi2v(i) end do close(u) ! initial beam data measured in mm -> transformed to cm x00v = 0.1_wp_*x00v y00v = 0.1_wp_*y00v z00v = 0.1_wp_*z00v waist1v = 0.1_wp_*waist1v waist2v = 0.1_wp_*waist2v rci1v = 10._wp_*rci1v rci2v = 10._wp_*rci2v iopt=0 call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier) call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier) call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier) call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier) call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier) call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier) call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier) call difcs(alphastv,x00v,nisteer,iopt,cx0,ier) call difcs(alphastv,y00v,nisteer,iopt,cy0,ier) call difcs(alphastv,z00v,nisteer,iopt,cz0,ier) if((alpha0 > alphastv(1)).and.(alpha0 < alphastv(nisteer))) then call locate(alphastv,nisteer,alpha0,k) dal=alpha0-alphastv(k) beta0=spli(cbeta,nisteer,k,dal) x00=spli(cx0,nisteer,k,dal) y00=spli(cy0,nisteer,k,dal) z00=spli(cz0,nisteer,k,dal) wcsi=spli(cwaist1,nisteer,k,dal) weta=spli(cwaist2,nisteer,k,dal) rcicsi=spli(crci1,nisteer,k,dal) rcieta=spli(crci2,nisteer,k,dal) phiw=spli(cphi1,nisteer,k,dal) phir=spli(cphi2,nisteer,k,dal) else ! alpha0 outside table range if(alpha0 >= alphastv(nisteer)) ii=nisteer if(alpha0 <= alphastv(1)) ii=1 beta0=betastv(ii) x00=x00v(ii) y00=y00v(ii) z00=z00v(ii) wcsi=waist1v(ii) weta=waist2v(ii) rcicsi=rci1v(ii) rcieta=rci2v(ii) phiw=phi1v(ii) phir=phi2v(ii) end if deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v, & phi1v,phi2v,x00v,y00v,z00v,cbeta, & cx0,cy0,cz0,cwaist1,cwaist2,crci1,crci2,cphi1,cphi2) end subroutine read_beam1 subroutine read_beam2(file_beam,beamid,alpha0,beta0,fghz,iox,x00,y00,z00, & wcsi,weta,rcicsi,rcieta,phiw,phir,unit) use utils, only : get_free_unit, intlin, locate use reflections, only : inside use dierckx, only : curfit, splev, surfit, bispev implicit none character(len=*), intent(in) :: file_beam integer, intent(in) :: beamid real(wp_), intent(inout) :: alpha0,beta0 real(wp_), intent(out) :: fghz,x00,y00,z00, wcsi,weta,rcicsi,rcieta,phir,phiw integer, intent(out) :: iox integer, intent(in), optional :: unit character(len=20) :: beamname integer :: u integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta integer :: iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk integer :: nxwaist1, nywaist1, nxwaist2, nywaist2, nxrci1, nyrci1, nxrci2 integer :: nyrci2, nxphi1, nyphi1, nxphi2, nyphi2, nxx0, nyx0, nxy0, nyy0 integer :: nxz0, nyz0, kx, ky, ii, npolyg, nmax, lwrk2, in integer :: nxycoord integer, DIMENSION(:), ALLOCATABLE :: iwrk real(wp_) :: alphast,betast, waist1, waist2, rci1, rci2, phi1, phi2 real(wp_) :: fp, minx, maxx, miny, maxy, eps, xcoord0, ycoord0 real(wp_), DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv, & betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord, & ycoord, wrk, txwaist1, tywaist1, txwaist2, tywaist2, & txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2, & txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1, & cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2, & xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC, & ypolygC, xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC real(wp_), DIMENSION(4) :: xvert, yvert real(wp_), dimension(1) :: fi integer, parameter :: kspl=1 real(wp_), parameter :: sspl=0.01_wp_ if (present(unit)) then u=unit else u = get_free_unit() end if open(unit=u,file=file_beam,status='OLD',action='READ') !======================================================================================= ! # of beams read(u,*) nbeam ! ! unused beams' data jumprow=0 ! c==================================================================================== do i=1,beamid-1 read(u,*) beamname, iox, fghz, nalpha, nbeta jumprow = jumprow+nalpha*nbeta end do ! c==================================================================================== ! ! beam of interest read(u,*) beamname, iox, fghz, nalpha, nbeta ! ! c==================================================================================== ! unused beams' data grids do i=1,(nbeam - beamid) read(u,*) beamname end do do i=1,jumprow read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 end do ! c==================================================================================== ! ! # of elements in beam data grid nisteer = nalpha*nbeta ! allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), & waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), & phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer), & xcoord(nisteer),ycoord(nisteer)) ! ! c==================================================================================== ! beam data grid reading do i=1,nisteer read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 ! ! initial beam data (x00, y00, z00) are measured in mm -> transformed to cm x00v(i)=0.1d0*x00 y00v(i)=0.1d0*y00 z00v(i)=0.1d0*z00 alphastv(i)=alphast betastv(i)=betast waist1v(i)=0.1d0*waist1 rci1v(i)=1.0d1*rci1 waist2v(i)=0.1d0*waist2 rci2v(i)=1.0d1*rci2 phi1v(i)=phi1 phi2v(i)=phi2 end do close(u) ! c==================================================================================== ! ! fdeg = 0 alpha, beta free variables ! 1 alpha free variable ! 2 beta free variable ! 3 no free variables fdeg = 2*(1/nalpha) + 1/nbeta !####################################################################################### ! ! no free variables if(fdeg.eq.3) then alpha0=alphastv(1) beta0=betastv(1) x00=x00v(1) y00=y00v(1) z00=z00v(1) wcsi=waist1v(1) weta=waist2v(1) rcicsi=rci1v(1) rcieta=rci2v(1) phiw=phi1v(1) phir=phi2v(1) return end if !####################################################################################### ! ! !####################################################################################### if(fdeg.eq.2) then ! beta = independent variable ! alpha = dependent variable xcoord = betastv ycoord = alphastv xcoord0 = beta0 ycoord0 = alpha0 kx=min(nbeta-1,kspl) ! c==================================================================================== else ! c==================================================================================== ! alpha = independent variable ! beta = dependent/independent (fdeg = 1/0) xcoord = alphastv ycoord = betastv xcoord0 = alpha0 ycoord0 = beta0 nxcoord = nalpha nycoord = nbeta kx=min(nalpha-1,kspl) ky=min(nbeta-1,kspl) end if !####################################################################################### ! iopt = 0 incheck = 0 ! !####################################################################################### if(fdeg.ne.0) then nxest = kx + nxcoord + 1 lwrk = (nxcoord*(kx+1)+nxest*(7+3*kx)) kwrk = nxest allocate(cycoord(nxest), txycoord(nxest), cwaist1(nxest), & txwaist1(nxest), cwaist2(nxest), txwaist2(nxest), & crci1(nxest), txrci1(nxest), crci2(nxest), txrci2(nxest), & cphi1(nxest), txphi1(nxest), cphi2(nxest), txphi2(nxest), & cx0(nxest), txx0(nxest), cy0(nxest), txy0(nxest), & cz0(nxest), txz0(nxest), w(nxcoord), wrk(lwrk), iwrk(kwrk)) ! w = 1.d0 ! ! 2D interpolation call curfit(iopt,nxcoord,xcoord,ycoord,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxycoord, & txycoord,cycoord,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,waist1v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist1, & txwaist1,cwaist1,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,waist2v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist2, & txwaist2,cwaist2,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,rci1v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci1, & txrci1,crci1,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,rci2v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci2, & txrci2,crci2,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,phi1v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi1, & txphi1,cphi1,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,phi2v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi2, & txphi2,cphi2,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,x00v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxx0, & txx0,cx0,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,y00v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxy0, & txy0,cy0,fp,wrk,lwrk,iwrk,ier) ! call curfit(iopt,nxcoord,xcoord,z00v,w, & xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxz0, & txz0,cz0,fp,wrk,lwrk,iwrk,ier) ! ! check if xcoord0 is out of table range ! incheck = 1 inside / 0 outside if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer)) incheck=1 ! c==================================================================================== else ! c==================================================================================== npolyg = 2*(nxcoord+nycoord-2) minx = minval(xcoord) maxx = maxval(xcoord) miny = minval(ycoord) maxy = maxval(ycoord) nxest = kx + 1 + int(sqrt(nisteer/2.)) nyest = ky + 1 + int(sqrt(nisteer/2.)) nmax = max(nxest,nyest) eps = 10.**(-8) lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19 lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2 kwrk = nisteer+(nmax-3)*(nmax-3) allocate(cwaist1(nxest*nyest), txwaist1(nmax), tywaist1(nmax), & cwaist2(nxest*nyest), txwaist2(nmax), tywaist2(nmax), & crci1(nxest*nyest), txrci1(nmax), tyrci1(nmax), & crci2(nxest*nyest), txrci2(nmax), tyrci2(nmax), & cphi1(nxest*nyest), txphi1(nmax), typhi1(nmax), & cphi2(nxest*nyest), txphi2(nmax), typhi2(nmax), & cx0(nxest*nyest), txx0(nmax), tyx0(nmax), & cy0(nxest*nyest), txy0(nmax), tyy0(nmax), & cz0(nxest*nyest), txz0(nmax), tyz0(nmax), & wrk(lwrk), wrk2(lwrk2), iwrk(kwrk), & xpolyg(npolyg), ypolyg(npolyg), w(nisteer)) ! w = 1.d0 ! ! 3D interpolation call surfit(iopt,nisteer,xcoord,ycoord,waist1v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxwaist1,txwaist1,nywaist1,tywaist1,cwaist1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,waist2v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxwaist2,txwaist2,nywaist2,tywaist2,cwaist2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,rci1v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxrci1,txrci1,nyrci1,tyrci1,crci1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,rci2v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxrci2,txrci2,nyrci2,tyrci2,crci2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,phi1v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxphi1,txphi1,nyphi1,typhi1,cphi1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,phi2v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxphi2,txphi2,nyphi2,typhi2,cphi2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,x00v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxx0,txx0,nyx0,tyx0,cx0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,y00v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxy0,txy0,nyy0,tyy0,cy0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! call surfit(iopt,nisteer,xcoord,ycoord,z00v,w, & minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxz0,txz0,nyz0,tyz0,cz0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! data range polygon xpolyg(1:nxcoord) = xcoord(1:nxcoord) ypolyg(1:nxcoord) = ycoord(1:nxcoord) ! ! c==================================================================================== do i=1,nycoord-2 xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord) xpolyg(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1) ypolyg(nxcoord+i) = ycoord((i+1)*nxcoord) ypolyg(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1) end do ! c==================================================================================== do i=1,nxcoord xpolyg(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1) ypolyg(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1) end do ! c==================================================================================== ! ! check if (xcoord0, ycoord0) is out of table range ! incheck = 1 inside / 0 outside if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0)) incheck = 1 end if deallocate(wrk,iwrk) !####################################################################################### ! ! !####################################################################################### if(fdeg.ne.0) then ! c==================================================================================== if(incheck.eq.1) then call splev(txycoord,nxycoord,cycoord,kx,(/xcoord0/),fi,1,ier) ycoord0=fi(1) call splev(txwaist1,nxwaist1,cwaist1,kx,(/xcoord0/),fi,1,ier) wcsi=fi(1) call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier) weta=fi(1) call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier) rcicsi=fi(1) call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier) rcieta=fi(1) call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier) phiw=fi(1) call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier) phir=fi(1) call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier) x00=fi(1) call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier) y00=fi(1) call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier) z00=fi(1) ! c---------------------------------------------------------------------------------- else ! c---------------------------------------------------------------------------------- if(xcoord0.ge.xcoord(nisteer)) ii=nisteer if(xcoord0.le.xcoord(1)) ii=1 ! xcoord0=xcoord(ii) ycoord0=ycoord(ii) x00=x00v(ii) y00=y00v(ii) z00=z00v(ii) wcsi=waist1v(ii) weta=waist2v(ii) rcicsi=rci1v(ii) rcieta=rci2v(ii) phiw=phi1v(ii) phir=phi2v(ii) end if ! c==================================================================================== else ! c==================================================================================== if(incheck.eq.0) then allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), & ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), & xpolygD(nycoord), ypolygD(nycoord), & xoutA(nxcoord+3), youtA(nxcoord+3), xoutB(nycoord+3), & youtB(nycoord+3), xoutC(nxcoord+3), youtC(nxcoord+3)) ! coordinates of vertices v1,v2,v3,v4 xvert(1) = xpolyg(1) xvert(2) = xpolyg(nxcoord) xvert(3) = xpolyg(nxcoord+nycoord-1) xvert(4) = xpolyg(2*nxcoord+nycoord-2) yvert(1) = ypolyg(1) yvert(2) = ypolyg(nxcoord) yvert(3) = ypolyg(nxcoord+nycoord-1) yvert(4) = ypolyg(2*nxcoord+nycoord-2) ! coordinates of side A,B,C,D xpolygA = xpolyg(1:nxcoord) ypolygA = ypolyg(1:nxcoord) xpolygB = xpolyg(nxcoord:nxcoord+nycoord-1) ypolygB = ypolyg(nxcoord:nxcoord+nycoord-1) xpolygC = xpolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) ypolygC = ypolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) xpolygD(1:nycoord-1) = xpolyg(2*nxcoord+nycoord-2:npolyg) xpolygD(nycoord) = xpolyg(1) ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg) ypolygD(nycoord) = ypolyg(1) ! contour of outside regions A (1,2), B(3,4), C(5,6) xoutA = huge(one) xoutA(1:nxcoord) = xpolygA xoutA(nxcoord+3) = xvert(1) youtA = -huge(one) youtA(1:nxcoord) = ypolygA youtA(nxcoord+1) = yvert(2) xoutB = huge(one) xoutB(1:nycoord) = xpolygB xoutB(nycoord+1) = xvert(3) youtB = huge(one) youtB(1:nycoord) = ypolygB youtB(nycoord+3) = yvert(2) xoutC = -huge(one) xoutC(1:nxcoord) = xpolygC xoutC(nxcoord+3) = xvert(3) youtC = huge(one) youtC(1:nxcoord) = ypolygC youtC(nxcoord+1) = yvert(4) ! c---------------------------------------------------------------------------------- ! search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid ! ! (6) | (5) | (4) ! _ _ _v4__________________v3 _ _ _ _ _ ! \ C \ ! \ \ (1)->(8) outside regions ! (7) D \ \ B (3) v1->v4 grid vertices ! \ \ A-D grid sides ! _ _ _ _ _ _\_________________\_ _ _ _ ! v1 A v2 ! (8) | (1) | (2) ! if(inside(xoutA,youtA,nxcoord+3,xcoord0,ycoord0)) then in = 1 if(xcoord0.GT.xvert(2)) then in = 2 end if else if(inside(xoutB,youtB,nycoord+3,xcoord0,ycoord0)) then in = 3 if(ycoord0.GT.yvert(3)) then in = 4 end if else if(inside(xoutC,youtC,nxcoord+3,xcoord0,ycoord0)) then in = 5 if(xcoord0.LT.xvert(4)) then in = 6 end if else in = 7 if(ycoord0.LT.yvert(1)) then in = 8 end if end if ! c---------------------------------------------------------------------------------- ! (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border ! depending on the region ! 1: xcoord0 unchanged, ycoord0 moved on side A ! 3: xcoord0 moved on side B, ycoord0 unchanged ! 5: xcoord0 unchanged, ycoord0 moved on side C ! 7: xcoord0 moved on side D, ycoord0 unchanged ! 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates ! in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi in ! new (xcoord0,ycoord0) ! in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the ! (xcoord0,ycoord0) vertex are used alpha0 = xcoord0 beta0 = ycoord0 SELECT CASE (in) CASE (1) ! beta0 outside table range ! locate position of xcoord0 with respect to x coordinates of side A call locate(xpolygA,nxcoord,xcoord0,ii) ! find corresponding y value on side A for xcoord position call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0) incheck = 1 CASE (2) ! alpha0 and beta0 outside table range ! xcoord0, ycoord0 set xcoord0 = xvert(2) ycoord0 = yvert(2) ii = nxcoord !indice per assegnare valori waist, rci, phi CASE (3) ! alpha0 outside table range call locate(ypolygB,nycoord,ycoord0,ii) call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0) incheck = 1 CASE (4) ! alpha0 and beta0 outside table range xcoord0 = xvert(3) ycoord0 = yvert(3) ii = nxcoord+nycoord-1 CASE (5) ! beta0 outside table range call locate(xpolygC,nxcoord,xcoord0,ii) call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0) incheck = 1 CASE (6) ! alpha0 and beta0 outside table range xcoord0 = xvert(4) ycoord0 = yvert(4) ii = 2*nxcoord+nycoord-2 CASE (7) ! alpha0 outside table range call locate(ypolygD,nycoord,ycoord0,ii) call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0) incheck = 1 CASE (8) ! alpha0 and beta0 outside table range xcoord0 = xvert(1) ycoord0 = yvert(1) ii = 1 END SELECT ! c---------------------------------------------------------------------------------- ! deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, xpolygD, ypolygD) end if ! c==================================================================================== ! ! c==================================================================================== if(incheck.eq.1) then lwrk = 2*(kx+ky+2) kwrk = 4 allocate(wrk(lwrk),iwrk(kwrk)) call bispev(txwaist1,nxwaist1,tywaist1,nywaist1,cwaist1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) wcsi=fi(1) call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) weta=fi(1) call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) rcicsi=fi(1) call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) rcieta=fi(1) call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) phiw=fi(1) call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) phir=fi(1) call bispev(txx0,nxx0,tyx0,nyx0,cx0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) x00=fi(1) call bispev(txy0,nxy0,tyy0,nyy0,cy0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) y00=fi(1) call bispev(txz0,nxz0,tyz0,nyz0,cz0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) z00=fi(1) deallocate(wrk,iwrk) ! c---------------------------------------------------------------------------------- else ! c---------------------------------------------------------------------------------- x00=x00v(ii) y00=y00v(ii) z00=z00v(ii) wcsi=waist1v(ii) weta=waist2v(ii) rcicsi=rci1v(ii) rcieta=rci2v(ii) phiw=phi1v(ii) phir=phi2v(ii) end if ! c==================================================================================== end if !####################################################################################### ! if(fdeg.ne.0) then deallocate(cycoord, txycoord, cwaist1, txwaist1, cwaist2, & txwaist2, crci1, txrci1, crci2, txrci2, cphi1, txphi1, & cphi2, txphi2, cx0, txx0, cy0, txy0, cz0, txz0, w) else deallocate(cwaist1, txwaist1, tywaist1, cwaist2, txwaist2, tywaist2, & crci1, txrci1, tyrci1, crci2, txrci2, tyrci2, & cphi1, txphi1, typhi1, cphi2, txphi2, typhi2, & cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0, & wrk2, xpolyg, ypolyg, w) end if ! !####################################################################################### ! set correct values for alpha, beta if(fdeg.eq.2) then alpha0 = ycoord0 beta0 = xcoord0 else alpha0 = xcoord0 beta0 = ycoord0 end if !####################################################################################### deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & phi2v,x00v,y00v,z00v,xcoord,ycoord) ! end subroutine read_beam2 subroutine launchangles2n(alpha,beta,xv,anv) use const_and_precisions, only : degree implicit none ! arguments real(wp_), intent(in) :: alpha,beta,xv(3) real(wp_), intent(out) :: anv(3) ! local variables real(wp_) :: r,anr,anphi,a,b r=sqrt(xv(1)**2+xv(2)**2) ! phi=atan2(y,x) a = degree*alpha b = degree*beta ! ! angles alpha, beta in a local reference system as proposed by Gribov et al ! anr = -cos(b)*cos(a) anphi = sin(b) ! anx = -cos(b)*cos(a) ! any = sin(b) anv(1) = (anr*xv(1) - anphi*xv(2))/r ! = anx anv(2) = (anr*xv(2) + anphi*xv(1))/r ! = any ! anr = (anx*xv(1) + any*xv(2))/r ! anphi = (any*xv(1) - anx*xv(2))/r anv(3) =-cos(b)*sin(a) ! = anz end subroutine launchangles2n subroutine xgygcoeff(fghz,ak0,bres,xgcn) use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_,pi,wce1_ implicit none ! arguments real(wp_), intent(in) :: fghz real(wp_), intent(out) :: ak0,bres,xgcn ! local variables real(wp_) :: omega omega=2.0e9_wp_*pi*fghz ! [rad/s] ak0=omega/vc ! [rad/cm] ! ! yg=btot/bres ! bres=omega/wce1_ ! [T] ! ! xg=xgcn*dens19 ! xgcn=4.0e13_wp_*pi*qe**2/(me*omega**2) ! [10^-19 m^3] end subroutine xgygcoeff end module beams