module beams use const_and_precisions, only : wp_, one implicit none contains subroutine read_beam0(params, err) ! Reads the wave launcher parameters for the simple case ! where w(z) and 1/R(z) are fixed. ! ! The file should be formatted as follows: ! ! 1 f ! 2 x₀ y₀ z₀ ! 3 w₀₁ w₀₂ d₀₁ d₀₂ φ ! ! where: ! - f is the frequency (GHz) ! - x₀, y₀, z₀ are the launcher position (cm) ! - w₀₁,w₀₂ are the beam waists in the two principal directions (cm) ! - d₀₁,d₀₂ are the distances of the beam waists from the launch ! point (cm) ! - φ is the rotation angle from the horizontal direction to the ! first principal direction (deg) ! ! Note: this case implies simple astigmatism, i.e. the ! amplitude and phase ellipses in the transverse plane ! are aligned. use const_and_precisions, only : pi, vc=>ccgs_ use gray_params, only : antenna_parameters use logger, only : log_error ! subroutine arguments type(antenna_parameters), intent(inout) :: params integer, intent(out) :: err ! local variables integer :: u real(wp_) :: k0, w0(2), d0(2), z_R(2), phi open(newunit=u, file=params%filenm, status='old', action='read', iostat=err) if (err /= 0) then call log_error('opening beams file ('//trim(params%filenm)//') failed!', & mod='beams', proc="read_beam0") return end if read(u, *) params%fghz ! Wave frequency (GHz) read(u, *) params%pos ! Launcher position (Cartesian coordinates) read(u, *) w0, & ! Beam waists (in the ξ,η axes) d0, & ! Waists (along ξ,η axes) distance from launcher phi ! rotation angle from horizontal direction to ! ξ axis in the transverse plane close(u) ! Wavevector k₀ k0 = 2.0e9_wp_* pi * params%fghz / vc ! Rayleigh ranges in the ξ, η directions: ! z_R = ½k₀w₀² z_R = k0/2 * w0**2 ! Beam widths in the ξ, η directions at z=0: ! w² = w₀² (1 + z₀²/z_R²) params%w = w0 * sqrt(1 + (d0/z_R)**2) ! Curvature in the ξ, η directions at z=0: ! K = 1/R = (z-z₀)/[(z-z₀)² + z_R²] params%ri = -d0 / (d0**2 + z_R**2) ! Ellipse axes angle params%phi = phi end subroutine read_beam0 subroutine read_beam1(params, err) ! Reads the wave launcher parameters for the case ! where w(z, α) and 1/R(z, α) depend on the launcher angle α. ! ! Format notes: ! 1. The first two lines contain the frequency in GHz ! and the number of rows. ! 2. The rest of file is a table with the following columns: ! θ, α, β, x₀, y₀, z₀, w₁, w₂, k₁, k₂, φ_w, φ_R ! 3 The meaning of the columns is ! - θ is the mechanical steering angle (unused) ! - α, β are the poloidal and toroidal launch angles (deg) ! - x₀, y₀, z₀ are the launcher position (mm) ! - w₁,w₂ are the beam widths in the two principal directions (mm) ! - k₁,k₂ are the wavefront curvatures in the two principal ! directions (mm⁻¹) ! - φ_w, φ_R are the rotation angles of the amplitude and phase ! ellipses in the transverse plane at the launch point (deg) use gray_params, only : antenna_parameters use splines, only : spline_simple use utils, only : locate use logger, only : log_error ! subroutine arguments type(antenna_parameters), intent(inout) :: params integer, intent(out) :: err ! local variables integer :: u, nisteer, i, k, ii real(wp_) :: steer, dal real(wp_), dimension(:), allocatable :: & alphastv, betastv, x00v, y00v, & z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v type(spline_simple) :: beta, waist1, waist2, & rci1, rci2, phi1, phi2, & x0, y0, z0 open(newunit=u, file=params%filenm, status='old', action='read', iostat=err) if (err /= 0) then call log_error('opening beams file ('//trim(params%filenm)//') failed!', & mod='beams', proc="read_beam1") return end if read(u,*) params%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)) 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 * rci1v rci2v = 10 * rci2v call beta%init(alphastv, betastv, nisteer) call waist1%init(alphastv, waist1v, nisteer) call waist2%init(alphastv, waist2v, nisteer) call rci1%init(alphastv, rci1v, nisteer) call rci2%init(alphastv, rci2v, nisteer) call phi1%init(alphastv, phi1v, nisteer) call phi2%init(alphastv, phi2v, nisteer) call x0%init(alphastv, x00v, nisteer) call y0%init(alphastv, y00v, nisteer) call z0%init(alphastv, z00v, nisteer) if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then call locate(alphastv, params%alpha, k) dal = params%alpha - alphastv(k) params%beta = beta%raw_eval(k, dal) params%pos(1) = x0%raw_eval(k, dal) params%pos(2) = y0%raw_eval(k, dal) params%pos(3) = z0%raw_eval(k, dal) params%w(1) = waist1%raw_eval(k, dal) params%w(2) = waist2%raw_eval(k, dal) params%ri(1) = rci1%raw_eval(k, dal) params%ri(2) = rci2%raw_eval(k, dal) params%phi(1) = phi1%raw_eval(k, dal) params%phi(2) = phi2%raw_eval(k, dal) else ! params%alpha outside table range if(params%alpha >= alphastv(nisteer)) ii=nisteer if(params%alpha <= alphastv(1)) ii=1 params%alpha = alphastv(ii) params%beta = betastv(ii) params%pos(1) = x00v(ii) params%pos(2) = y00v(ii) params%pos(3) = z00v(ii) params%w(1) = waist1v(ii) params%w(2) = waist2v(ii) params%ri(1) = rci1v(ii) params%ri(2) = rci2v(ii) params%phi(1) = phi1v(ii) params%phi(2) = phi2v(ii) end if deallocate(alphastv, betastv, waist1v, waist2v, & rci1v, rci2v, phi1v, phi2v, & x00v, y00v, z00v) call beta%deinit call waist1%deinit call waist2%deinit call rci1%deinit call rci2%deinit call phi1%deinit call phi2%deinit call x0%deinit call y0%deinit call z0%deinit end subroutine read_beam1 subroutine read_beam2(params, beamid, err) ! Reads the wave launcher parameters for the general case ! where w(z, α, β) and 1/R(z, α, β) depend on the launcher angles α, β. ! ! Format notes: ! 1. The first line specifies the number N of beams described by the file. ! 2. The rest of the files consists of N 2D tables preceded by a header ! id, mode, f, na, nb ! where ! - id is a string identifier of the beam ! - mode indicates where the beam has O-mode (1) or X-mode (2) ! polarisation ! - f is the frequency (GHz) ! - nα, nβ are the numbers of rows and columns of the table. ! 3. The 2D table is stored in row-major order over *nα×nβ* lines, that is, the ! i,j-th record is stored on the l-th line, with l = i + nα*j. ! 4. The poloidal angle *α(i,j)* must be monotonic along *i* and the toroidal angle ! *β(i, j)* must be monotonic along *j*. ! 5. Each line stores one record with the same fields as in the 1D format. use gray_params, only : antenna_parameters use utils, only : linear_interp, locate use types, only : contour use dierckx, only : curfit, splev, surfit, bispev use logger, only : log_error ! subroutine arguments type(antenna_parameters), intent(inout) :: params integer, intent(in) :: beamid integer, intent(out), optional :: err ! local variables 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 type(contour) :: polyg, polygA, polygB, polygC, polygD, outA, outB, outC real(wp_), DIMENSION(4) :: xvert, yvert real(wp_), dimension(1) :: fi integer, parameter :: kspl=1 real(wp_), parameter :: sspl=0.01_wp_ open(newunit=u, file=params%filenm, status='old', action='read', iostat=err) if (err /= 0) then call log_error('opening beams file ('//trim(params%filenm)//') failed!', & mod='beams', proc="read_beam1") return end if !======================================================================================= ! # of beams read(u,*) nbeam ! ! unused beams' data jumprow=0 ! c==================================================================================== do i=1,beamid-1 read(u,*) beamname, params%iox, params%fghz, nalpha, nbeta jumprow = jumprow+nalpha*nbeta end do ! c==================================================================================== ! ! beam of interest read(u,*) beamname, params%iox, params%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,params%pos(1),params%pos(2),params%pos(3),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,params%pos(1),params%pos(2),params%pos(3),waist1,waist2,rci1,rci2,phi1,phi2 ! ! initial beam data (params%pos(1), params%pos(2), params%pos(3)) are measured in mm -> transformed to cm x00v(i)=0.1d0*params%pos(1) y00v(i)=0.1d0*params%pos(2) z00v(i)=0.1d0*params%pos(3) 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 params%alpha=alphastv(1) params%beta=betastv(1) params%pos(1)=x00v(1) params%pos(2)=y00v(1) params%pos(3)=z00v(1) params%w(1)=waist1v(1) params%w(2)=waist2v(1) params%ri(1)=rci1v(1) params%ri(2)=rci2v(1) params%phi(2)=phi1v(1) params%phi(1)=phi2v(1) deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & phi2v,x00v,y00v,z00v,xcoord,ycoord) return end if !####################################################################################### ! ! !####################################################################################### if(fdeg.eq.2) then ! beta = independent variable ! alpha = dependent variable xcoord = betastv ycoord = alphastv xcoord0 = params%beta ycoord0 = params%alpha kx=min(nbeta-1,kspl) ! c==================================================================================== else ! c==================================================================================== ! alpha = independent variable ! beta = dependent/independent (fdeg = 1/0) xcoord = alphastv ycoord = betastv xcoord0 = params%alpha ycoord0 = params%beta 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 = 2*(kx + 1) nyest = 2*(ky + 1) nmax = max(nxest,nyest)+max(kx,ky) 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), & polyg%R(npolyg), polyg%z(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 polyg%R(1:nxcoord) = xcoord(1:nxcoord) polyg%z(1:nxcoord) = ycoord(1:nxcoord) ! ! c==================================================================================== do i=1,nycoord-2 polyg%R(nxcoord+i) = xcoord((i+1)*nxcoord) polyg%R(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1) polyg%z(nxcoord+i) = ycoord((i+1)*nxcoord) polyg%z(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1) end do ! c==================================================================================== do i=1,nxcoord polyg%R(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1) polyg%z(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 (polyg%contains(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) params%w(1)=fi(1) call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier) params%w(2)=fi(1) call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier) params%ri(1)=fi(1) call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier) params%ri(2)=fi(1) call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier) params%phi(2)=fi(1) call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier) params%phi(1)=fi(1) call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier) params%pos(1)=fi(1) call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier) params%pos(2)=fi(1) call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier) params%pos(3)=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) params%pos(1)=x00v(ii) params%pos(2)=y00v(ii) params%pos(3)=z00v(ii) params%w(1)=waist1v(ii) params%w(2)=waist2v(ii) params%ri(1)=rci1v(ii) params%ri(2)=rci2v(ii) params%phi(2)=phi1v(ii) params%phi(1)=phi2v(ii) end if ! c==================================================================================== else ! c==================================================================================== if(incheck.eq.0) then allocate(polygA%R(nxcoord), polygA%z(nxcoord)) allocate(polygC%R(nxcoord), polygC%z(nxcoord)) allocate(polygB%R(nycoord), polygB%z(nycoord)) allocate(polygD%R(nycoord), polygD%z(nycoord)) allocate(outA%R(nxcoord+3), outA%z(nxcoord+3)) allocate(outB%R(nycoord+3), outB%z(nycoord+3)) allocate(outC%R(nxcoord+3), outC%z(nxcoord+3)) ! coordinates of vertices v1,v2,v3,v4 xvert(1) = polyg%R(1) xvert(2) = polyg%R(nxcoord) xvert(3) = polyg%R(nxcoord+nycoord-1) xvert(4) = polyg%R(2*nxcoord+nycoord-2) yvert(1) = polyg%z(1) yvert(2) = polyg%z(nxcoord) yvert(3) = polyg%z(nxcoord+nycoord-1) yvert(4) = polyg%z(2*nxcoord+nycoord-2) ! coordinates of side A,B,C,D polygA%R = polyg%R(1:nxcoord) polygA%z = polyg%z(1:nxcoord) polygB%R = polyg%R(nxcoord:nxcoord+nycoord-1) polygB%z = polyg%z(nxcoord:nxcoord+nycoord-1) polygC%R = polyg%R(nxcoord+nycoord-1:2*nxcoord+nycoord-2) polygC%z = polyg%z(nxcoord+nycoord-1:2*nxcoord+nycoord-2) polygD%R(1:nycoord-1) = polyg%R(2*nxcoord+nycoord-2:npolyg) polygD%R(nycoord) = polyg%R(1) polygD%z(1:nycoord-1) = polyg%z(2*nxcoord+nycoord-2:npolyg) polygD%z(nycoord) = polyg%z(1) ! contour of outside regions A (1,2), B(3,4), C(5,6) outA%R = huge(one) outA%R(1:nxcoord) = polygA%R outA%R(nxcoord+3) = xvert(1) outA%z = -huge(one) outA%z(1:nxcoord) = polygA%z outA%z(nxcoord+1) = yvert(2) outB%R = huge(one) outB%R(1:nycoord) = polygB%R outB%R(nycoord+1) = xvert(3) outB%z = huge(one) outB%z(1:nycoord) = polygB%z outB%z(nycoord+3) = yvert(2) outC%R = -huge(one) outC%R(1:nxcoord) = polygC%R outC%R(nxcoord+3) = xvert(3) outC%z = huge(one) outC%z(1:nxcoord) = polygC%z outC%z(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 (outA%contains(xcoord0, ycoord0)) then in = 1 if(xcoord0.GT.xvert(2)) then in = 2 end if else if (outB%contains(xcoord0, ycoord0)) then in = 3 if(ycoord0.GT.yvert(3)) then in = 4 end if else if (outC%contains(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 params%pos(1),params%pos(2),params%pos(3),waist,rci,phi in ! new (xcoord0,ycoord0) ! in 2,4,6,8 incheck remains 0 and params%pos(1),params%pos(2),params%pos(3),waist,rci,phi values at the ! (xcoord0,ycoord0) vertex are used params%alpha = xcoord0 params%beta = ycoord0 SELECT CASE (in) CASE (1) ! params%beta outside table range ! locate position of xcoord0 with respect to x coordinates of side A call locate(polygA%R, xcoord0, ii) ! find corresponding y value on side A for xcoord position ycoord0 = linear_interp(polygA%R(ii), polygA%z(ii), polygA%R(ii+1), polygA%z(ii+1), xcoord0) incheck = 1 CASE (2) ! params%alpha and params%beta outside table range ! xcoord0, ycoord0 set xcoord0 = xvert(2) ycoord0 = yvert(2) ii = nxcoord !indice per assegnare valori waist, rci, phi CASE (3) ! params%alpha outside table range call locate(polygB%z, ycoord0, ii) xcoord0 = linear_interp(polygB%z(ii), polygB%R(ii), polygB%z(ii+1), polygB%R(ii+1), ycoord0) incheck = 1 CASE (4) ! params%alpha and params%beta outside table range xcoord0 = xvert(3) ycoord0 = yvert(3) ii = nxcoord+nycoord-1 CASE (5) ! params%beta outside table range call locate(polygC%R, xcoord0, ii) ycoord0 = linear_interp(polygC%R(ii+1), polygC%z(ii+1), polygC%R(ii), polygC%z(ii), xcoord0) incheck = 1 CASE (6) ! params%alpha and params%beta outside table range xcoord0 = xvert(4) ycoord0 = yvert(4) ii = 2*nxcoord+nycoord-2 CASE (7) ! params%alpha outside table range call locate(polygD%z, ycoord0, ii) xcoord0 = linear_interp(polygD%z(ii), polygD%R(ii), polygD%z(ii+1), polygD%R(ii+1), ycoord0) incheck = 1 CASE (8) ! params%alpha and params%beta outside table range xcoord0 = xvert(1) ycoord0 = yvert(1) ii = 1 END SELECT ! c---------------------------------------------------------------------------------- ! deallocate(polygA%R, polygA%z, polygC%R, polygC%z, polygB%R, polygB%z, & polygD%R, polygD%z, outA%R, outA%z, outB%R, outB%z, outC%R, outC%z) 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) params%w(1)=fi(1) call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%w(2)=fi(1) call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%ri(1)=fi(1) call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%ri(2)=fi(1) call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%phi(2)=fi(1) call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%phi(1)=fi(1) call bispev(txx0,nxx0,tyx0,nyx0,cx0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%pos(1)=fi(1) call bispev(txy0,nxy0,tyy0,nyy0,cy0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%pos(2)=fi(1) call bispev(txz0,nxz0,tyz0,nyz0,cz0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) params%pos(3)=fi(1) deallocate(wrk,iwrk) ! c---------------------------------------------------------------------------------- else ! c---------------------------------------------------------------------------------- params%pos(1)=x00v(ii) params%pos(2)=y00v(ii) params%pos(3)=z00v(ii) params%w(1)=waist1v(ii) params%w(2)=waist2v(ii) params%ri(1)=rci1v(ii) params%ri(2)=rci2v(ii) params%phi(2)=phi1v(ii) params%phi(1)=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, polyg%R, polyg%z, w) end if ! !####################################################################################### ! set correct values for alpha, beta if(fdeg.eq.2) then params%alpha = ycoord0 params%beta = xcoord0 else params%alpha = xcoord0 params%beta = ycoord0 end if !####################################################################################### deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & phi2v,x00v,y00v,z00v,xcoord,ycoord) ! end subroutine read_beam2 subroutine launchangles2n(params, N) ! Given the wave launcher `params` computes the initial ! wavevector N̅ = ck̅/ω. ! ! Notes: ! 1. the result is in Cartesian coordinates; ! 2. the beam is assumed to start in a vacuum, so N = 1. use const_and_precisions, only : degree use gray_params, only : antenna_parameters ! subroutine arguments type(antenna_parameters), intent(in) :: params real(wp_), intent(out) :: N(3) ! local variables real(wp_) :: r, Nr, Nphi, a, b R = hypot(params%pos(1), params%pos(2)) a = params%alpha * degree b = params%beta * degree ! N̅ in cylindrical coordinates using the poloidal ! and toroidal launch angles. Nr = -cos(a)*cos(b) Nphi = sin(b) ! Convert to Cartesian coordinates ! Notes: ! 1. The unit vectors in Cartesian coordinates are: ! R^ = x^ cosφ + y^ sinφ ! φ^ = x^ sinφ - y^ sinφ ! 2. tanφ = y/x ⇒ { cosφ = x/√(x²+y²) = x/R ! { sinφ = y/√(x²+y²) = x/R N(1) = (Nr*params%pos(1) - Nphi*params%pos(2)) / R N(2) = (Nr*params%pos(2) + Nphi*params%pos(1)) / R N(3) = -cos(b)*sin(a) end subroutine launchangles2n subroutine xgygcoeff(fghz, k0, Bres, xgcn) ! Given the EC wave frequency computes: ! ! 1. vacuum wavevector `k0` (k₀ = ω/c), ! 2. resonant magnetic field `Bres` (qB/m = ω), ! 3. adimensional `xgcn` parameter (X = ω_p²/ω² = nq²/ε₀mω²). use const_and_precisions, only : qe=>ecgs_, me=>mecgs_, & vc=>ccgs_, pi, wce1_ ! subroutine arguments real(wp_), intent(in) :: fghz real(wp_), intent(out) :: k0, Bres, xgcn ! local variables real(wp_) :: omega omega = 2.0e9_wp_*pi*fghz ! [rad/s] k0 = 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 pure subroutine gaussian_eikonal(Q, r, z, S, grad_S, hess_S) ! Computes the complex eikonal, gradient and Hessian of a Gaussian beam ! subroutine arguments complex(wp_), intent(in) :: Q(2,2) ! complex curvature tensor Q(z) real(wp_), intent(in) :: r(2), z ! coordinates complex(wp_), intent(out) :: S ! complex eikonal complex(wp_), intent(out) :: grad_S(3) ! gradient of imaginary eikonal complex(wp_), intent(out) :: hess_S(3,3) ! Hessian of imaginary eikonal ! local variables complex(wp_) :: Q2(2,2), Q3(2,2) Q2 = matmul(Q, Q) Q3 = matmul(Q, Q2) ! The most general Gaussian beam is given by ! ! E̅(r̅,z) = e̅(z) exp[-ik₀z - ik₀r̅⋅Q(z)⋅r̅/2 + iη(z)] ! ! where: ! - e̅(z) is the amplitude ! - Q(z) = Q₀ (I + zQ₀)⁻¹ is the complex curvature tensor ! - η(z) is the Gouy phase ! - k₀ = ω/c ! - r̅ = [x, y] ! ! See https://doi.org/10.1364/AO.52.006030 for more details. ! ! Since the eikonal ansatz in GRAY is E̅(x̅,t) = e̅(r̅) exp[-ik₀S(x̅) + iωt], ! the eikonal for a Gaussian beam is: ! ! S(x̅) = z + ½ r̅⋅Q(z)⋅r̅ - η(z)/k₀ ! S = z + dot_product(r, matmul(Q, r))/2 ! Splitting up x̅ into (r̅, z), the gradient of S can be written as ! ! ∇S(r̅,z) = [Q(z)⋅r̅, 1 + ½ r̅⋅dQ/dz⋅r̅] ! ! For dQ/dz, we use the matrix identity dA/dz = - A⁻¹ dA/dz A⁻¹: ! ! dQ/dz = Q₀ d/dz (I + zQ₀)⁻¹ ! = - Q₀ (I + zQ₀)⁻¹ Q₀ (I + zQ₀)⁻¹ ! = - Q(z)² ! ! So ∇S(r̅,z) = [Q(z)⋅r̅, 1 - ½ r̅⋅Q(z)²⋅r̅] ! grad_S = [matmul(Q, r), 1 - dot_product(r, matmul(Q2, r))/2] ! Carrying on, the Hessian splits into 4 blocks (2x2, 2x1, 1x2, 1x1): ! ! HS(r̅,z) = k₀ [ Q(z), dQ/dz⋅r̅] = k₀ [ Q(z), -Q(z)²⋅r̅] ! [-Q(z)²⋅r̅, -½ r̅⋅dQ²/dz⋅r̅] [-Q(z)²⋅r̅, r̅⋅Q³⋅r̅] ! ! having used dQ²/dz = 2Q dQ/dz = -2Q³. ! hess_S(1:2,:) = reshape([Q, -matmul(Q2, r)], [2,3]) hess_S(3, :) = [-matmul(Q2, r), dot_product(r, matmul(Q3, r))] ! Note: this ignores the Gouy phase, but for the computation of the ! initial conditions (x̅, N̅_R=∇S_R, N̅_I=∇S_I for each ray) it is ! inconsequential because: ! ! 1. the initial conditions are given at z=const, so η=const as well ! 2. N̅_I = ∇S_I does not depend on η(z) ! 3. N̅_R does depend on η(z), but is obtained indirectly as the ! solution of the vacuum quasioptical dispersion relation: ! ! { N̅_R ⋅ N̅_I = 0 ! { N_R² - N_I² = N₀² ! ! The latter implies N̅_I and N̅_R are not exactly consistent with ! a Gaussian beam, but this is a necessary trade-off. end subroutine gaussian_eikonal end module beams