subroutine gray_integration use beamdata, only : nstep implicit none integer istep,istop,index_rt real*8 st,dst,strfl11 integer i real*8 st0 common/ss/st common/dsds/dst common/istep/istep common/istop/istop common/strfl11/strfl11 common/index_rt/index_rt c ray integration: begin st0=0.0d0 if(index_rt.gt.1) st0=strfl11 do i=1,nstep istep=i st=i*dst+st0 c advance one step call rkint4 c calculations after one step: call after_onestep(istep,istop) if(istop.eq.1) exit c end do c ray integration: end return end subroutine after_gray_integration(rhopin,nrho,dpdvout,ajcdout) use beamdata, only : nrayr implicit real*8 (a-h,o-z) parameter(zero=0.0d0) character*255 filenmeqq,filenmprf,filenmbm dimension rhopin(nrho),dpdvout(nrho),ajcdout(nrho) c common/ss/st common/ibeam/ibeam common/warm/iwarm,ilarm common/filesn/filenmeqq,filenmprf,filenmbm common/iieq/iequil common/iipr/iprof common/index_rt/index_rt c common/p0/p0mw common/factb/factb common/taumnx/taumn,taumx,pabstot,currtot common/scal/iscal common/facttn/factt,factn common/pardens/dens0,aln1,aln2 common/parqte/te0,dte0,alt1,alt2 c c print all ray positions in local reference system c if(nrayr.gt.1) then iproj=1 nfilp=609 call projxyzt(iproj,nfilp) end if c c print final results on screen c write(*,*) write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st if(iwarm.gt.0) then write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx else write(*,'(a,2f9.4)') 'taumn, taumx = ', zero,zero end if c write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot currtka =currtot*1.0d3 write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka c if (index_rt.eq.1) then if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq if(iequil.eq.1) write(*,*) 'ANALTYCAL EQUILIBRIUM' if(iprof.eq.1) write(*,*) 'PROFILE file : ',filenmprf if(iprof.eq.0) write(*,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0 if(ibeam.ge.1) write(*,*) 'LAUNCHER CASE : ',filenmbm end if c c compute power and current density profiles for all rays c pabs=pabstot currt=currtot call pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout) c return 99 format(20(1x,e16.8e3)) end c c c subroutine after_onestep(i,istop) use const_and_precisions, only : pi use beamdata, only : psjki,ppabs,ccci,iiv,tauv, . iop,iow,tau1v,yyrfl,nrayr,nrayth implicit real*8 (a-h,o-z) parameter(taucr=12.0d0,cvdr=pi/180.0d0) dimension xv(3),anv(3),xvrfl(3),anvrfl(3) common/warm/iwarm,ilarm common/ist/istpr0,istpl0 common/istgr/istpr,istpl c common/psinv/psinv common/psinv11/psinv11 common/ierr/ierr common/taumnx/taumn,taumx,pabstot,currtot common/xv/xv common/anv/anv common/cent/btrcen,rcen c common/p0/p0mw common/pol0/psipol0,chipol0 common/ipol/ipolc common/iovmin/iopmin,iowmin common/densbnd/psdbnd common/powrfl/powrfl common/dstvac/dstvac common/strfl11/strfl11 common/dsds/dst common/index_rt/index_rt common/ipass/ipass common/rwallm/rwallm common/zbound/zbmin,zbmax pabstot=0.0d0 currtot=0.0d0 taumn=1d+30 taumx=-1d+30 psinv11=1.0d0 iopmin=100 iowmin=100 c do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk call gwork(j,k) c if(ierr.gt.0) then write(*,*) ' IERR = ', ierr if(ierr.eq.97) then c igrad=0 c ierr=0 else istop=1 c exit end if end if psjki(j,k,i)=psinv rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0 zzm=xv(3)/100.d0 if(j.eq.1) rrm11=rrm if (iwarm.gt.0.and.i.gt.1) then if(psinv.ge.0.and.psinv.le.1.0d0.and. . zzm.ge.zbmin.and.zzm.le.zbmax) then call pabs_curr(i,j,k) iiv(j,k)=i else if(iiv(j,k).gt.1) tauv(j,k,i)=tauv(j,k,i-1) end if if(tauv(j,k,i).le.taumn) taumn=tauv(j,k,i) if(tauv(j,k,i).ge.taumx) taumx=tauv(j,k,i) pabstot=pabstot+ppabs(j,k,iiv(j,k)) currtot=currtot+ccci(j,k,iiv(j,k)) end if call print_output(i,j,k) if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) . iop(j,k)=1 if(iop(j,k).eq.1.and. . (psinv.ge.psdbnd.or. . (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax)))) . iop(j,k)=2 c iov=0 initially, iov=1 first entrance in plasma, c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma if(index_rt.eq.1) then if(j.eq.1) then psinv11=psinv if(ipolc.eq.0.and.iop(j,k).eq.1) then call pol_limit(qqin,uuin,vvin) ipolc=1 qqa=cos(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr) uua=sin(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr) vva=sin(2.0d0*chipol0*cvdr) powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin) c p0mw=p0mw*powa c write(*,*) ' ' c write(*,*) 'Coupled power fraction =',powa c write(*,*) ' ' c write(*,*) 'Input coupled power (MW) =',p0mw c write(*,*) ' ' end if if (iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1 . .and.ipolc.eq.1) then call pol_limit(qqout,uuout,vvout) ipolc=2 call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) strfl11=i*dst+dstvac call pol_limit(qqin2,uuin2,vvin2) if(irfl.gt.0) then powrfl=0.5d0*(1.0d0+vvrfl*vvin2+uurfl*uuin2+qqrfl*qqin2) else powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2) end if write(*,*) 'Reflected power fraction =',powrfl iop(j,k)=3 yyrfl(j,k,1)=xvrfl(1) yyrfl(j,k,2)=xvrfl(2) yyrfl(j,k,3)=xvrfl(3) yyrfl(j,k,4)=anvrfl(1) yyrfl(j,k,5)=anvrfl(2) yyrfl(j,k,6)=anvrfl(3) tau1v(j,k)=tauv(j,k,iiv(j,k)) end if else if(iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1) then call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) iop(j,k)=3 yyrfl(j,k,1)=xvrfl(1) yyrfl(j,k,2)=xvrfl(2) yyrfl(j,k,3)=xvrfl(3) yyrfl(j,k,4)=anvrfl(1) yyrfl(j,k,5)=anvrfl(2) yyrfl(j,k,6)=anvrfl(3) tau1v(j,k)=tauv(j,k,iiv(j,k)) end if end if end if if(iop(j,k).lt.iopmin) iopmin=iop(j,k) end do end do if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1 psimin=psjki(1,1,i) if(nrayr.gt.1) . psimin=min(psimin,minval(psjki(2:nrayr,1:nrayth,i))) if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1) . istop=1 if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1 if(iopmin.eq.2.and.ipass.eq.1) istop=1 if(iopmin.eq.3) istop=1 c print ray positions for j=nrayr in local reference system istpr=istpr+1 if (istpr.eq.istpr0) then c write(*,*) 'istep = ',i if(nrayr.gt.1) then iproj=0 nfilp=608 call projxyzt(iproj,nfilp) end if istpr=0 end if c if (istpl.eq.istpl0) istpl=0 istpl=istpl+1 return end c c c subroutine print_output(i,j,k) use const_and_precisions, only : pi use beamdata, only : ywrk,psjki,tauv,alphav,pdjki, . currj,didst,q,tau1v,nrayr!,nrayth implicit real*8 (a-h,o-z) parameter(taucr=12.0d0) complex*16 ex,ey,ez c common/nhn/nhn common/iokh/iohkawa common/p0/p0mw common/index_rt/index_rt common/ss/st common/istgr/istpr,istpl common/ist/istpr0,istpl0 common/iieq/iequil c common/parpl/brr,bphi,bzz,ajphi common/btot/btot common/xgxg/xg common/ygyg/yg common/dddens/dens,ddens common/tete/tekev common/absor/alpha,effjcd,akim,tau0 common/densbnd/psdbnd common/epolar/ex,ey,ez c common/nplr/anpl,anpr common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nprw/anprr,anpri c x=ywrk(1,j,k) y=ywrk(2,j,k) z=ywrk(3,j,k) rr=sqrt(x*x+y*y) c anx=ywrk(4,j,k) any=ywrk(5,j,k) anz=ywrk(6,j,k) anr=(anx*x+any*y)/rr anphi=(any*x-anx*y)/rr c cnphi=(any*x-anx*y) c rrm=rr*1.0d-2 zzm=z*1.0d-2 stm=st*1.0d-2 xxm=x*1.0d-2 yym=y*1.0d-2 c if(index_rt.gt.1) then taujk=tauv(j,k,i)+tau1v(j,k) else taujk=tauv(j,k,i) end if c central ray only begin if(j.eq.1) then phi=acos(x/sqrt(y*y+x*x)) if(y.lt.0.0d0) phi=-phi phideg=phi*180.0d0/pi psi=psjki(j,k,i) rhot=1.0d0 bbr=0.0d0 bbz=0.0d0 bpol=0.0d0 rhot=1.0d0 dens11=0.0d0 if(psi.ge.0.0d0) then if(iequil.eq.2) then if (psi.le.1.0d0) rhot=frhotor(psi) bbr=brr bbz=bzz bpol=sqrt(bbr**2+bbz**2) else rhot=sqrt(psi) end if else tekev=0.0d0 akim=0.0d0 end if ddr11=ddr c cutoff=xg-(1-yg)*(1-anpl**2) c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 if(psi.le.psdbnd.and.psi.ge.0.0d0) dens11=dens dids11=didst(j,k,i)*1.0d2/(p0mw*q(j)) dpdv11=pdjki(j,k,i)/(p0mw*q(j)) ajcd11=currj(j,k,i)/(p0mw*q(j)) alpha11=alphav(j,k,i) tau11=taujk pt11=exp(-taujk) write(604,99) stm,rrm,zzm,phideg, . psi,rhot,dens11,tekev, . bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100, . alpha11,tau11,pt11,ajcd11,dids11, . dble(nhn),dble(iohkawa),dble(index_rt) c call polarcold(exf,eyif,ezf,elf,etf) end if c central ray only end if(k.eq.1.and.j.eq.nrayr) write(617,99) st,ddr11,ddr,ddi c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3 if(istpl.eq.istpl0) then if(j.eq.nrayr) then c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then write(633,111) i,j,k,stm,xxm,yym,rrm,zzm, . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) end if c if(k.eq.nrayth) write(633,*) ' ' end if c return 99 format(30(1x,e16.8e3)) 111 format(3i5,16(1x,e16.8e3)) end c c c subroutine prfile implicit none integer index_rt common/index_rt/index_rt If(index_rt.eq.1) then write(604,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '// .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' write(608,*) ' #istep j k xt yt zt rt psin' write(609,*) ' #istep j k xt yt zt rt psin' write(617,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(633,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' write(612,*) ' #i sst psi w1 w2' write(607,*)'#Icd Pa Jphip dPdVp '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' write(648,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' else c If(index_rt.eq.3) then write(604,*) ' ' write(608,*) ' ' write(609,*) ' ' write(617,*) ' ' write(612,*) ' ' write(648,*) ' ' c end if end if return end c c c subroutine read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, . dne, zeff, qsf, powin, alphain, betain, beamid) use green_func_p, only:Setup_SpitzFunc use const_and_precisions, only : pi use beamdata, only : nrayr,nrayth,nstep implicit real*8 (a-h,o-z) integer ijetto, mr, mz, nrho, nbnd, beamid real*8 r(mr), z(mz), psin(mr,mz) real*8 psiax, psibnd, rax, zax, powin real*8 alphain, betain real*8 rbnd(nbnd), zbnd(nbnd) real*8 psijet(nrho), f(nrho), te(nrho), dne(nrho), . zeff(nrho), qsf(nrho) real*8 me character*8 wdat character*10 wtim character*255 filenmeqq,filenmprf,filenmbm parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(cvdr=pi/180.0d0) parameter(nbb=5000) real*8 rlim(nbb),zlim(nbb) c common/xgcn/xgcn common/ipec/ipec,nnd common/ibeam/ibeam common/ist/istpr0,istpl0 common/warm/iwarm,ilarm common/ieccd/ieccd common/idst/idst common/imx/imx c common/filesn/filenmeqq,filenmprf,filenmbm c common/rwmax/rwmax common/dsds/dst common/igrad/igrad common/ipass/ipass common/rwallm/rwallm common/limiter/rlim,zlim,nlim common/iieq/iequil common/icocos/icocos common/ixp/ixp common/ipsn/ipsinorm common/sspl/sspl common/iipr/iprof common/factb/factb common/facttn/factt,factn common/cent/btrcen,rcen c common/parwv/ak0,akinv,fhz common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/pol0/psipol0,chipol0 c common/pardens/dens0,aln1,aln2 common/parban/b0,rr0m,zr0m,rpam common/parqq/q0,qa,alq common/parqte/te0,dte0,alt1,alt2 common/zz/Zeffan c common/parbres/bres common/densbnd/psdbnd common/nfile/neqdsk,nprof common/sgnib/sgnbphi,sgniphi common/p0/p0mw c common/mode/sox common/angles/alpha0,beta0 common/scal/iscal common/fghz/fghz c open(602,file='gray.data',status= 'unknown') c c alpha0, beta0 (cartesian) launching angles c fghz wave frequency (GHz) c p0mw injected power (MW) c read(602,*) dummy, dummy !alpha0,beta0 c angles read from values passed as arguments (in rad) alpha0 = alphain beta0 = betain read(602,*) fghz read(602,*) p0mw c power value overwritten with value passed as argument (in W) p0mw=powin*1.d-6 c c nrayr number of rays in radial direction c nrayth number of rays in angular direction c rwmax normalized maximum radius of beam power c rwmax=1 -> last ray at radius = waist c read(602,*) nrayr,nrayth,rwmax if(nrayr.eq.1) nrayth=1 c c x00,y00,z00 coordinates of launching point c read(602,*) x00,y00,z00 c c beams parameters in local reference system c w0 -> waist, d0 -> waist distance from launching point c phiw angle of beam ellipse c read(602,*) w0csi,w0eta,d0csi,d0eta,phiw c c ibeam=0 :read data for beam as above c ibeam=1 :read data from file simple astigmatic beam c ibeam=2 :read data from file general astigmatic beam read(602,*) ibeam read(602,*) filenmbm c c iequil=0 :vacuum c iequil=1 :analytical equilibrium c iequil=2 :read eqdsk c ixp=0,-1,+1 : no X point , bottom/up X point c read(602,*) iequil,ixp c c iprof=0 :analytical density and temp. profiles c iprof>0 :numerical density and temp. profiles c read(602,*) iprof c c iwarm=0 :no absorption and cd c iwarm=1 :weakly relativistic absorption c iwarm=2 :relativistic absorption, n<1 asymptotic expansion c iwarm=3 :relativistic absorption, numerical integration c ilarm :order of larmor expansion c imx :max n of iterations in dispersion, imx<0 uses 1st c iteration in case of failure after |imx| iterations c read(602,*) iwarm,ilarm,imx c if(iwarm.gt.2) iwarm=2 c c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models c read(602,*) ieccd if (ieccd.ge.10) call Setup_SpitzFunc c c ipec=0/1/-1 :pec profiles grid equispaced in psi/rhop/as input c nnd :number of grid steps for pec profiles +1 c read(602,*) ipec,nnd c input radial grid forced on output. nnd is ignored in this case ipec=-1 c c dst integration step c nstep maximum number of integration steps c istpr0 projection step = dsdt*istprj c istpl0 plot step = dsdt*istpl c igrad=0 optical ray-tracing, initial conditions as for beam c igrad=1 quasi-optical ray-tracing c igrad=-1 ray-tracing, init. condit. c from center of mirror and with angular spread c ipass=1/2 1 or 2 passes into plasma c iox=1/2 OM/XM c idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr c read(602,*) dst,nstep,istpr0,istpl0,idst read(602,*) igrad,ipass,rwallm read(602,*) iox,psipol0,chipol0 c c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c sspl spline parameter for psi interpolation c read(602,*) filenmeqq read(602,*) ipsinorm,sspl c psin grid is normalized ipsinorm=1 c c factb factor for magnetic field (only for numerical equil) c scaling adopted: beta=const, qpsi=const, nustar=const c iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling c factT factn factor for Te&ne scaling c read(602,*) factb, iscal,factt,factn if (factb.eq.0.0d0) factb=1.0d0 c c psbnd value of psi ( > 1 ) of density boundary c read(602,*) filenmprf read(602,*) psdbnd if(iprof.eq.0) psdbnd=1.0d0 c c signum of toroidal B and I c icocos index for equilibrium from COCOS - O. Sauter Feb 2012 read(602,*) sgnbphi,sgniphi,icocos c data are already in the coordinate system used by gray icocos=3 c c read parameters for analytical density and temperature c profiles if iprof=0 c if(iprof.eq.0) then read(602,*) dens0,aln1,aln2 read(602,*) te0,dte0,alt1,alt2 else read(602,*) dummy,dummy,dummy read(602,*) dummy,dummy,dummy,dummy end if read(602,*) zeffan c c read parameters for analytical simple cylindical equilibrium c if iequil=0,1 if(iequil.lt.2) then read(602,*) rr0,zr0,rpa read(602,*) b0 read(602,*) q0,qa,alq rr0m=rr0/1.0d2 zr0m=zr0/1.0d2 rpam=rpa/1.0d2 else read(602,*) dummy,dummy,dummy read(602,*) dummy read(602,*) dummy,dummy,dummy end if c close(unit=602) c if(nrayr.eq.1) igrad=0 if (nrayr.lt.5) then igrad=0 write(*,*) ' nrayr < 5 ! => OPTICAL CASE ONLY' write(*,*) end if c fhz=fghz*1.0d9 ak0=2.0d9*pi*fghz/vc akinv=1.0d0/ak0 c bresg=2.0d0*pi*fhz*me*vc/qe bres=bresg*1.0d-4 c c xg=xgcn*dens19 c xgcn=1.0d-5*qe**2/(pi*me*fghz**2) c c read data for beam from file if ibeam>0 c if(ibeam.gt.0) then filenmbm='graybeam.data' call read_beams(beamid,iox) else zrcsi=0.5d0*ak0*w0csi**2 zreta=0.5d0*ak0*w0eta**2 if(igrad.gt.0) then wcsi=w0csi*sqrt(1.0d0+(d0csi/zrcsi)**2) weta=w0eta*sqrt(1.0d0+(d0eta/zreta)**2) rcicsi=-d0csi/(d0csi**2+zrcsi**2) rcieta=-d0eta/(d0eta**2+zreta**2) phir=phiw else d0eta=d0csi wcsi=w0csi*abs(d0csi/zrcsi) weta=w0eta*abs(d0eta/zreta) rcicsi=w0csi/zrcsi rcieta=w0eta/zreta if(d0csi.gt.0) then rcicsi=-rcicsi rcieta=-rcieta end if phir=phiw end if end if c sox=-1.0d0 if(iox.eq.2) sox=1.0d0 c write(*,'(a,2f8.3)') 'alpha0, beta0 = ',alpha0,beta0 c r00=sqrt(x00**2+y00**2) phi0=acos(x00/r00) if(y00.lt.0) phi0=-phi0 write(*,'(a,4f8.3)') 'x00, y00, R00, z00 = ',x00,y00,r00,z00 write(*,*) c c anx0c=(anr0c*x00-anphi0c*y00)/r00 c any0c=(anr0c*y00+anphi0c*x00)/r00 c c anr0c=(anx0c*x00+any0c*y00)/r00 c anphi0c=(any0c*x00-anx0c*y00)/r00 c c angles alpha0, beta0 in a local reference system as proposed by Gribov et al c c anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0) c anphi0c=sin(cvdr*beta0) c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0) anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0) anphi0c=sin(cvdr*beta0) anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0) anx0c=(anr0c*x00-anphi0c*y00)/r00 any0c=(anr0c*y00+anphi0c*x00)/r00 c c read data for Te , ne , Zeff from file if iprof>0 c if (iprof.eq.1) then c nprof=98 c lprf=length(filenmprf) c open(file=filenmprf(1:lprf)//'.prf', c . status= 'unknown',unit=nprof) c read profiles from input arguments filenmprf='JETTO' call profdata(nrho, psijet, te, dne, zeff) c close(nprof) end if c c read equilibrium data from file if iequil=2 c if (iequil.eq.2) then c neqdsk=99 c leqq=length(filenmeqq) c open(file=filenmeqq(1:leqq)//'.eqdsk', c . status= 'unknown', unit=neqdsk) filenmeqq='JETTO' call equidata(ijetto, mr, mz, r, z, psin, psiax, psibnd, . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, qsf) c close(neqdsk) c print density, temperature, safecty factor, toroidal current dens c versus psi, rhop, rhot call print_prof end if if (iequil.eq.1) call surf_anal if (iequil.ne.2.or.ipass.lt.0) then c set simple limiter as two cylindrical walls at rwallm and r00 nlim=5 rlim(1)=rwallm rlim(2)=max(r00*1.d-2,rmxm) rlim(3)=rlim(2) rlim(4)=rlim(1) rlim(5)=rlim(1) zlim(1)=(z00-dst*nstep*2)*1.d-2 zlim(2)=zlim(1) zlim(3)=(z00+dst*nstep*2)*1.d-2 zlim(4)=zlim(3) zlim(5)=zlim(1) ipass=abs(ipass) end if nfil=638 open(nfil,file='headers.txt',status= 'unknown') call date_and_time(wdat,wtim) write(nfil,916) wdat(1:4),wdat(5:6),wdat(7:8), . wtim(1:2),wtim(3:4),wtim(5:6) write(nfil,904) REVISION if (iequil.eq.2) then write(nfil,905) trim(filenmeqq) else if (iequil.eq.1) write(nfil,912) rr0,zr0,rpa,B0,q0,qa,alq if (iequil.eq.0) write(nfil,'("# VACUUM CASE")') end if if (iprof.eq.1) then write(nfil,907) trim(filenmprf) else write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeffan end if write(nfil,911) fghz,p0mw,iox write(nfil,902) x00,y00,z00 write(nfil,908) alpha0,beta0 if(ibeam.ge.1) write(nfil,909) trim(filenmbm) if(ibeam.eq.0) write(nfil,903) w0csi,w0eta,d0csi,d0eta,phiw write(nfil,900) nrayr, nrayth, rwmax write(nfil,901) igrad,iwarm,ilarm,ieccd,idst write(nfil,915) dst,nstep write(nfil,914) sgnbphi,sgniphi,icocos write(nfil,906) factb,factt,factn,iscal write(nfil,910) sspl,psdbnd,nnd,ipec,ipsinorm write(nfil,*) ' ' close(nfil) return 900 format('# Nray_r Nray_th rwmax : ',2i5,1x,es12.5) 901 format('# igrad iwarm ilarm ieccd idst : ',6i5) 902 format('# X0 Y0 Z0 launching point (cm) : ',3(1x,es12.5)) 903 format('# w0xi w0eta d0xi d0eta (cm) phiw (deg) : ',5(1x,es12.5)) 904 format('# GRAY revision : ',a) 905 format('# EQUILIBRIUM file : ',a) 906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5)) 907 format('# PROFILES file : ',a) 908 format('# alpha0 beta0 launch angles (deg) CYL : ',2(1x,es12.5)) 909 format('# LAUNCHER file : ',a24) 910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5) 911 format('# fghz P0 IOX : ',2(1x,es12.5),i5) 912 format('# AN. EQUILIBRIUM R0 z0 a B0 q0 qa alq : ',7(1x,es12.5)) 913 format('# AN. PROFILES ne0 aln1 aln2 Te0 Tea alt1 alt2 Zeff : ' . ,8(1x,es12.5)) 914 format('# sgnB_phi sgnI_phi icocos : ',2(1x,es12.5),i5) 915 format('# dst nstep : ',1x,es12.5,i5) 916 format('# Date : ',a4,2('/',a2),1x,a2,2(':',a2)) 99 format(20(1x,e16.8e3)) end c c c subroutine surf_anal use const_and_precisions, only : pi implicit real*8(a-h,o-z) common/parban/b0,rr0m,zr0m,rpam common/parbres/bres c c print circular magnetic surfaces iequil=1 c write(631,*) '#i psin R z' dal=2.0d-2*pi drn=0.1d0 do i=1,10 rni=i*drn do k=1,101 drrm=rpam*rni*cos((k-1)*dal) dzzm=rpam*rni*sin((k-1)*dal) rrm=rr0m+drrm zzm=zr0m+dzzm write(631,111) i,rni,rrm,zzm end do write(631,*) ' ' write(631,*) ' ' end do c c print resonant B iequil=1 c write(630,*)'#i Btot R z' rres=b0*rr0m/bres zmx=zr0m+rpam zmn=zr0m-rpam do i=1,21 zzres=zmn+(i-1)*(zmx-zmn)/2.0d1 write(630,111) i,bres,rres,zzres end do return 111 format(i4,20(1x,e16.8e3)) end c c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ subroutine read_beams(beamid,iox) use reflections, only : inside implicit real*8(a-h,o-z) c character*255 filenmeqq,filenmprf,filenmbm character*20 beamname c integer beamid, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta, . iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk, . nxwaist1, nywaist1, nxwaist2, nywaist2, nxrci1, nyrci1, nxrci2, . nyrci2, nxphi1, nyphi1, nxphi2, nyphi2, nxx0, nyx0, nxy0, nyy0, . nxz0, nyz0, kx, ky, ii, nfbeam, iox, npolyg, nmax, lwrk2, in integer, DIMENSION(:), ALLOCATABLE :: iwrk c real*8 alphast, betast, x00, y00, z00, waist1, waist2, dvir1, . dvir2, delta, phi1, phi, zr1, zr2, ako, w1, w2, rci1, rci2, fghz, . fhz, wcso, weta, rcicsi, rcieta, phiw, phir, xcoord0, ycoord0, . alpha0, beta0, fp, minx, maxx, miny, maxy, eps, dmin real*8, DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv, . betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord, . ycoord, ycoord2, 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 real*8 xvert(4), yvert(4) c parameter(kspl=3,sspl=0.01) c common/ibeam/ibeam common/filesn/filenmeqq,filenmprf,filenmbm common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/mirr/x00,y00,z00 common/angles/alpha0,beta0 common/parwv/ak0,akinv,fhz common/fghz/fghz c c for given alpha0 -> beta0 + beam parameters c c ibeam=1 simple astigmatic beam c ibeam=2 general astigmatic beam c nfbeam=603 lbm=length(filenmbm) open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam) c c####################################################################################### if(ibeam.eq.1) then read(nfbeam,*) nisteer naplha = nisteer nbeta = 1 allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), . waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), . phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer)) c c c==================================================================================== do i=1,nisteer read(nfbeam,*) alphast,betast,x00,y00,z00, waist1,dvir1, . waist2,dvir2,delta phi1=delta phi2=delta zr1=0.5d-1*ak0*waist1**2 zr2=0.5d-1*ak0*waist2**2 w1=waist1*sqrt(1.0d0+(dvir1/zr1)**2) w2=waist2*sqrt(1.0d0+(dvir2/zr2)**2) rci1=-dvir1/(dvir1**2+zr1**2) rci2=-dvir2/(dvir2**2+zr2**2) c c 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*w1 rci1v(i)=1.0d1*rci1 waist2v(i)=0.1d0*w2 rci2v(i)=1.0d1*rci2 phi1v(i)=phi1 phi2v(i)=phi2 end do c c==================================================================================== c c free variable = alpha fdeg = 1 c======================================================================================= else c======================================================================================= c # of beams read(nfbeam,*) nbeam c c unused beams' data jumprow=0 c c==================================================================================== do i=1,beamid-1 read(nfbeam,*) beamname, iox, fghz, nalpha, nbeta jumprow = jumprow+nalpha*nbeta end do c c==================================================================================== c c beam of interest read(nfbeam,*) beamname, iox, fghz, nalpha, nbeta fhz=fghz*1.0d9 c c c==================================================================================== c unused beams' data grids do i=1,(jumprow + (nbeam - beamid)) read(nfbeam,*) end do c c==================================================================================== c c # of elements in beam data grid nisteer = nalpha*nbeta c allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), . waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), . phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer)) c c c==================================================================================== c beam data grid reading do i=1,nisteer read(nfbeam,*) alphast,betast,x00,y00,z00, . waist1,waist2,rci1,rci2,phi1,phi2 c c 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 c c==================================================================================== c c fdeg = 0 alpha, beta free variables c 1 alpha free variable c 2 beta free variable c 3 no free variables fdeg = 2*1/nalpha + 1/nbeta end if c####################################################################################### c close(nfbeam) c c####################################################################################### c 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 c####################################################################################### c c c####################################################################################### if(fdeg.eq.2) then c beta = independent variable c alpha = dependent variable xcoord = betastv ycoord = alphastv xcoord0 = beta0 ycoord0 = alpha0 kx=min(nbeta-1,kspl) c c==================================================================================== else c c==================================================================================== c alpha = independent variable c 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 c####################################################################################### c iopt = 0 incheck = 0 c c####################################################################################### 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)) c do i=1,nxcoord w(i) = 1.d0 end do c c 2D interpolation call dierckx_curfit(iopt,nxcoord,xcoord,ycoord,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxycoord, . txycoord,cycoord,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,waist1v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist1, . txwaist1,cwaist1,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,waist2v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist2, . txwaist2,cwaist2,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,rci1v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci1, . txrci1,crci1,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,rci2v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci21, . txrci2,crci2,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,phi1v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi1, . txphi1,cphi1,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,phi2v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi2, . txphi2,cphi2,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,x00v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxx0, . txx0,cx0,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,y00v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxy0, . txy0,cy0,fp,wrk,lwrk,iwrk,ier) c call dierckx_curfit(iopt,nxcoord,xcoord,z00v,w, . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxz0, . txz0,cz0,fp,wrk,lwrk,iwrk,ier) c c check if xcoord0 is out of table range c incheck = 1 inside / 0 outside if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer)) . incheck=1 c c==================================================================================== else c c==================================================================================== npolyg = 2*(nxcoord+nycoord-2) minx = minval(xcoord) maxx = maxval(xcoord) miny = minval(ycoord) maxy = maxval(ycoord) nxest = kx + 1 + sqrt(nisteer/2.) nyest = ky + 1 + 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)) c do i=1,nisteer w(i) = 1.d0 end do c c 3D interpolation call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c call dierckx_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,lwkr2,iwrk,kwrk,ier) c data range polygon xpolyg(1:nxcoord) = xcoord(1:nxcoord) ypolyg(1:nxcoord) = ycoord(1:nxcoord) c c 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 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 c==================================================================================== c c check if (xcoord0, ycoord0) is out of table range c incheck = 1 inside / 0 outside if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0)) . incheck = 1 deallocate(wrk,iwrk) end if c####################################################################################### c c c####################################################################################### if(fdeg.ne.0) then c c==================================================================================== if(incheck.eq.1) then call dierckx_splev(txycoord,nxycoord,cycoord,kx,xcoord0, . ycoord0,1,ier) c call dierckx_splev(txycoord,nxycoord,cwaist1,kx,xcoord0,wcsi, . 1,ier) c call dierckx_splev(txwaist2,nxwaist2,cwaist2,kx,xcoord0,weta, . 1,ier) c call dierckx_splev(txrci1,nxrci1,crci1,kx,xcoord0,rcicsi, . 1,ier) c call dierckx_splev(txrci2,nxrci2,crci2,kx,xcoord0,rcieta, . 1,ier) c call dierckx_splev(txphi1,nxphi1,cphi1,kx,xcoord0,phiw,1,ier) c call dierckx_splev(txphi2,nxphi2,cphi2,kx,xcoord0,phir,1,ier) c call dierckx_splev(txx0,nxx0,cx0,kx,xcoord0,x00,1,ier) c call dierckx_splev(txy0,nxy0,cy0,kx,xcoord0,y00,1,ier) c call dierckx_splev(txz0,nxz0,cz0,kx,xcoord0,z00,1,ier) c c---------------------------------------------------------------------------------- else c c---------------------------------------------------------------------------------- write(*,*) ' alpha0 or beta0 outside table range !!!' if(xcoord0.ge.xcoord(nisteer)) ii=nisteer if(xcoord0.le.xcoord(1)) ii=1 c 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 c==================================================================================== else c c==================================================================================== if(incheck.eq.0) then allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), . ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), . xpolygD(nycoord), ypolygD(nycoord)) c 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) c 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) c c---------------------------------------------------------------------------------- c search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid c c | | c (6) (5) (4) c | | c _ _ _ v4 _________________v3_ _ _ _ c | C | (1)->(8) outside regions c | | c (7) D | | B (3) v1->v4 grid vertices c | | c _ _ _ _ |_________________| _ _ _ _ A-D grid sides c v1 A v2 c | | c (8) (1) (2) c | | c if(xcoord0.gt.xvert(1).and.xcoord0.lt.xvert(2).and. . ycoord0.le.maxval(ypolygA)) then in=1 else if(ycoord0.gt.yvert(2).and.ycoord0.lt.yvert(3).and. . xcoord0.ge.minval(xpolygB)) then in=3 else if(xcoord0.lt.xvert(3).and.xcoord0.gt.xvert(4).and. . ycoord0.ge.minval(ypolygC)) then in=5 else if(ycoord0.lt.yvert(4).and.ycoord0.gt.yvert(1).and. . xcoord0.le.maxval(xpolygD)) then in=7 else if(xcoord0.ge.xvert(2).and.ycoord0.le.yvert(2)) then in=2 else if(xcoord0.ge.xvert(3).and.ycoord0.ge.yvert(3)) then in=4 else if(xcoord0.le.xvert(4).and.ycoord0.ge.yvert(4)) then in=6 else if(xcoord0.le.xvert(1).and.ycoord0.le.yvert(1)) then in=8 endif c c---------------------------------------------------------------------------------- c c c---------------------------------------------------------------------------------- c (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border c depending on the region c 1: xcoord0 unchanged, ycoord0 moved on side A c 3: xcoord0 moved on side B, ycoord0 unchanged c 5: xcoord0 unchanged, ycoord0 moved on side C c 7: xcoord0 moved on side D, ycoord0 unchanged c 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates c in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi in c new (xcoord0,ycoord0) c in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the c (xcoord0,ycoord0) vertex are used alpha0 = xcoord0 beta0 = ycoord0 SELECT CASE (in) CASE (1) write(*,*) ' beta0 outside table range !!!' c locate position of xcoord0 with respect to x coordinates of side A call vlocate(xpolygA,nxcoord,xcoord0,ii) c 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) write(*,*) ' alpha0 and beta0 outside table range !!!' c xcoord0, ycoord0 set xcoord0 = xvert(2) ycoord0 = yvert(2) ii = nxcoord !indice per assegnare valori waist, rci, phi CASE (3) write(*,*) ' alpha0 outside table range !!!' call vlocate(ypolygB,nycoord,ycoord0,ii) call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1), . xpolygB(ii+1),ycoord0,xcoord0) incheck = 1 CASE (4) write(*,*) ' alpha0 and beta0 outside table range !!!' xcoord0 = xvert(3) ycoord0 = yvert(3) ii = nxcoord+nycoord-1 CASE (5) write(*,*) ' beta0 outside table range !!!' call vlocate(xpolygC,nxcoord,xcoord0,ii) call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii), . ypolygC(ii),xcoord0,ycoord0) incheck = 1 CASE (6) write(*,*) ' alpha0 and beta0 outside table range !!!' xcoord0 = xvert(4) ycoord0 = yvert(4) ii = 2*nxcoord+nycoord-2 CASE (7) write(*,*) ' alpha0 outside table range !!!' call vlocate(ypolygD,nycoord,ycoord0,ii) call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1), . xpolygD(ii+1),ycoord0,xcoord0) incheck = 1 CASE (8) write(*,*) ' alpha0 and beta0 outside table range !!!!' xcoord0 = xvert(1) ycoord0 = yvert(1) ii = 1 END SELECT c c---------------------------------------------------------------------------------- c c print new (alpha0, beta0) write(*,*) ' (alpha0, beta0) values changed' write(*,10) ' old (alpha0,beta0) = (', alpha0, . ',', beta0, ')' write(*,10) ' new (alpha0,beta0) = (', xcoord0, . ',', ycoord0, ')' 10 format(a27,f10.5,a1,f10.5,a1) end if c c==================================================================================== c c c==================================================================================== if(incheck.eq.1) then lwrk = 2*(kx+ky+2) kwrk = 4 allocate(wrk(lwrk),iwrk(kwrk)) c call dierckx_bispev(txwaist1,nxwaist1,tywaist1,nywaist1, . cwaist1,kx,ky,(/xcoord0/),1,(/ycoord0/),1,wcsi, . wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txwaist2,nxwaist2,tywaist2,nywaist2, . cwaist2,kx,ky,(/xcoord0/),1,(/ycoord0/),1,weta, . wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1,kx, . ky,(/xcoord0/),1,(/ycoord0/),1,rcicsi,wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2,kx, . ky,(/xcoord0/),1,(/ycoord0/),1,rcieta,wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1,kx, . ky,(/xcoord0/),1,(/ycoord0/),1,phiw,wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2,kx, . ky,(/xcoord0/),1,(/ycoord0/),1,phir,wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txx0,nxx0,tyx0,nyx0,cx0,kx,ky, . (/xcoord0/),1,(/ycoord0/),1,x00,wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txy0,nxy0,tyy0,nyy0,cy0,kx,ky, . (/xcoord0/),1,(/ycoord0/),1,y00,wrk,lwrk,iwrk,kwrk,ier) c call dierckx_bispev(txz0,nxz0,tyz0,nyz0,cz0,kx,ky, . (/xcoord0/),1,(/ycoord0/),1,z00,wrk,lwrk,iwrk,kwrk,ier) c c---------------------------------------------------------------------------------- else c 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 c==================================================================================== end if c####################################################################################### c c c####################################################################################### c set correct values for alpha, beta if(fdeg.eq.2) then alpha0 = ycoord0 beta0 = xcoord0 else alpha0 = xcoord0 beta0 = ycoord0 end if c####################################################################################### c return end c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ subroutine read_beams_OLD C implicit real*8(a-h,o-z) C character*255 filenmeqq,filenmprf,filenmbm C parameter(nstrmx=50) Cc C dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4) C dimension x00v(nstrmx),y00v(nstrmx),z00v(nstrmx) C dimension cx0(nstrmx,4),cy0(nstrmx,4),cz0(nstrmx,4) C dimension waist1v(nstrmx),waist2v(nstrmx) C dimension rci1v(nstrmx),rci2v(nstrmx) C dimension cwaist1(nstrmx,4),cwaist2(nstrmx,4) C dimension crci1(nstrmx,4),crci2(nstrmx,4) C dimension phi1v(nstrmx),phi2v(nstrmx) C dimension cphi1(nstrmx,4),cphi2(nstrmx,4) Cc C common/ibeam/ibeam C common/filesn/filenmeqq,filenmprf,filenmbm C common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir C common/mirr/x00,y00,z00 C common/angles/alpha0,beta0 C common/parwv/ak0,akinv,fhz Cc Cc for given alpha0 -> beta0 + beam parameters Cc Cc ibeam=1 simple astigmatic beam Cc ibeam=2 general astigmatic beam Cc Cc initial beam data are measured in mm -> transformed to cm Cc C nfbeam=603 C lbm=length(filenmbm) C open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam) Cc C read(nfbeam,*) nisteer C do i=1,nisteer C if(ibeam.eq.1) then C read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, C . waist1,dvir1,waist2,dvir2,delta C phi1=delta C phi2=delta C zr1=0.5d-1*ak0*waist1**2 C zr2=0.5d-1*ak0*waist2**2 C w1=waist1*sqrt(1.0d0+(dvir1/zr1)**2) C w2=waist2*sqrt(1.0d0+(dvir2/zr2)**2) C rci1=-dvir1/(dvir1**2+zr1**2) C rci2=-dvir2/(dvir2**2+zr2**2) C else C read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, C . w1,w2,rci1,rci2,phi1,phi2 C end if C x00v(i)=0.1d0*x00mm C y00v(i)=0.1d0*y00mm C z00v(i)=0.1d0*z00mm C alphastv(i)=alphast C betastv(i)=betast C waist1v(i)=0.1d0*w1 C rci1v(i)=1.0d1*rci1 C waist2v(i)=0.1d0*w2 C rci2v(i)=1.0d1*rci2 C phi1v(i)=phi1 C phi2v(i)=phi2 C end do Cc C iopt=0 C call difcsg(alphastv,betastv,nisteer,iopt,cbeta,ier) C call difcsg(alphastv,waist1v,nisteer,iopt,cwaist1,ier) C call difcsg(alphastv,rci1v,nisteer,iopt,crci1,ier) C call difcsg(alphastv,waist2v,nisteer,iopt,cwaist2,ier) C call difcsg(alphastv,rci2v,nisteer,iopt,crci2,ier) C call difcsg(alphastv,phi1v,nisteer,iopt,cphi1,ier) C call difcsg(alphastv,phi2v,nisteer,iopt,cphi2,ier) C call difcsg(alphastv,x00v,nisteer,iopt,cx0,ier) C call difcsg(alphastv,y00v,nisteer,iopt,cy0,ier) C call difcsg(alphastv,z00v,nisteer,iopt,cz0,ier) Cc C if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then C call vlocate(alphastv,nisteer,alpha0,k) C dal=alpha0-alphastv(k) C betst=fspli(cbeta,nstrmx,k,dal) C x00=fspli(cx0,nstrmx,k,dal) C y00=fspli(cy0,nstrmx,k,dal) C z00=fspli(cz0,nstrmx,k,dal) C wcsi=fspli(cwaist1,nstrmx,k,dal) C weta=fspli(cwaist2,nstrmx,k,dal) C rcicsi=fspli(crci1,nstrmx,k,dal) C rcieta=fspli(crci2,nstrmx,k,dal) C phiw=fspli(cphi1,nstrmx,k,dal) C phir=fspli(cphi2,nstrmx,k,dal) C else C write(*,*) ' alpha0 outside table range !!!' C if(alpha0.ge.alphastv(nisteer)) ii=nisteer C if(alpha0.le.alphastv(1)) ii=1 C betst=betastv(ii) C x00=x00v(ii) C y00=y00v(ii) C z00=z00v(ii) C wcsi=waist1v(ii) C weta=waist2v(ii) C rcicsi=rci1v(ii) C rcieta=rci2v(ii) C phiw=phi1v(ii) C phir=phi2v(ii) C end if C beta0=betst Cc C close(nfbeam) Cc C return end c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ c subroutine equidata(ijetto,mr,mz,r,z,psin2d,psiaxis,psiedge, . rax,zax,nbnd,rbnd,zbnd,mrho,psijet,fpjet,qjet) use const_and_precisions, only : pi use reflections, only : inside implicit real*8 (a-h,o-z) integer ijetto,mr,mz,nbnd,mrho real*8 r(mr),z(mz),psin2d(mr,mz) real*8 psiaxis,psiedge,rax,zax real*8 rbnd(nbnd),zbnd(nbnd) real*8 psijet(mrho),fpjet(mrho),qjet(mrho) parameter(nnw=501,nnh=501) parameter(nbb=5000) dimension fpol(nnw),qpsi(nnw) dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw) dimension rlim(nbb),zlim(nbb) c parameter(nrest=nnw+4,nzest=nnh+4) parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54) parameter(liwrk=nnh+nnw+nrest+nzest+3,kspl=3) dimension fvpsi(nnw*nnh),cc(nnw*nnh),ffvpsi(nnw*nnh) dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) c parameter(lwrkbsp=4*(nnw+nnh),liwrkbsp=nnw+nnh) c dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) dimension rv1d(nnw*nnh),zv1d(nnw*nnh),wpsi(nnw*nnh) parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) parameter(lw11=nnw*3+nnh*3+nnw*nnh,ldiwrk=nnw+nnh) dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) dimension iwrkd(ldiwrk) parameter(lwrkf=nnw*4+nrest*16) dimension tfp(nrest),cfp(nrest),wrkf(lwrkf),iwrkf(nrest),wf(nnw) dimension fpoli(nnw) c common/pareq1/psia common/pareq1t/psiant,psinop common/cent/btrcen,rcen common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/psinr/psinr common/qpsi/qpsi common/cpsin/rv,zv,psin common/cpsi/psi common/eqnn/nr,nz,nrho,npp,nintp common/ipsn/ipsinorm common/sspl/sspl common/nfile/neqdsk,nprof common/zbound/zbmin,zbmax common/sgnib/sgnbphi,sgniphi common/factb/factb common/ixp/ixp common/icocos/icocos common/coffeqt/tr,tz common/coffeqtp/tfp common/coffeq/cc common/coffeqd/cc01,cc10,cc20,cc02,cc11 common/coffeqn/nsrt,nszt,nsft common/cofffp/cfp common/fpas/fpolas common/rhotsx/rhotsx common/phitedge/phitedge common/rrtor/rrtor common/cnt/rup,zup,rlw,zlw common/limiter/rlim,zlim,nlim c c read from file eqdsk c (see http://fusion.gat.com/efit/g_eqdsk.html) c c spline interpolation of psi and derivatives c nr=mr nz=mz nrho=mrho do i=1,nr rv(i)=r(i) enddo do j=1,nz zv(j)=z(j) do i=1,nr psin(i,j)=psin2d(i,j) enddo enddo do i=1,nrho psinr(i)=psijet(i) fpol(i)=fpjet(i) qpsi(i)=qjet(i) end do rmaxis=rax zmaxis=zax c assume last point of fpol array corresponds to vacuum value c and compute vacuum field at magnetic axis as reference rcen=rax btrcen=fpol(nrho)/rcen sgnbphi=sign(1.0d0,btrcen) c c length in m !!! c dr=rv(2)-rv(1) dz=zv(2)-zv(1) rmnm=rv(1) rmxm=rv(nr) zmnm=zv(1) zmxm=zv(nz) c psi function psia0=psiedge-psiaxis psia=psia0*factb sgniphi=sign(1.0d0,-psia0) c c do j=1,nz c do i=1,nr c write(620,2021) rv(i),zv(j),psin(i,j) c enddo c write(620,*) ' ' c enddo c c spline fitting/interpolation of psin(i,j) and derivatives c if (ijetto.eq.1) then c valid data on the whole grid do j=1,nz do i=1,nr c psi(i,j)=psin(i,j)*psia ij=nz*(i-1)+j fvpsi(ij)=psin(i,j) enddo enddo iopt=0 call dierckx_regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, . wrk,lwrk,iwrk,liwrk,ier) c if ier=-1 data are fitted using sspl=0 if(ier.eq.-1) then sspl=0.0d0 call dierckx_regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm, . zmxm,kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc, . fp,wrk,lwrk,iwrk,liwrk,ier) end if nsrt=nsr nszt=nsz else c valid data only inside boundary ij=0 do j=1,nz do i=1,nr if (.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle ij=ij+1 rv1d(ij)=rv(i) zv1d(ij)=zv(j) fvpsi(ij)=psin(i,j) wpsi(ij)=1.0d0 end do end do c add boundary to the fit c c *** TEMPORARILY COMMENTED TO TRY REDUCING MEMORY CONSUMPTION - BEGIN *** c if (ij+nbnd.le.nnw*nnh) then c do i=1,nbnd c ij=ij+1 c rv1d(ij)=rbnd(i) c zv1d(ij)=zbnd(i) c fvpsi(ij)=1.0d0 c wpsi(ij)=1.0d0 c end do c end if c *** TEMPORARILY COMMENTED TO TRY REDUCING MEMORY CONSUMPTION - END *** c nrz=ij c c fit as a scattered set of points c nsr=nr/4+4 nsz=nz/4+4 call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cc,ier) c if ier=-1 data are fitted using sspl=0 if(ier.eq.-1) then sspl=0.0d0 nsr=nr/4+4 nsz=nz/4+4 call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cc,ier) end if nsrt=nsr nszt=nsz end if c c re-evaluate psi on the grid using the spline (only for debug and cniteq) c c call dierckx_bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi, c . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier) c c do j=1,nz c do i=1,nr c ij=nz*(i-1)+j c psin(i,j)=ffvpsi(ij) cc psi(i,j)=psin(i,j)*psia c write(619,2021) rv(i),zv(j),psin(i,j) c enddo c write(619,*) ' ' c enddo c c2021 format(5(1x,e16.9)) nur=1 nuz=0 call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv, . nz,ffvpsi,cc10,lw10,iwrkd,ldiwrk,ier) c nur=0 nuz=1 call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv, . nz,ffvpsi,cc01,lw01,iwrkd,ldiwrk,ier) c nur=2 nuz=0 call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv, . nz,ffvpsi,cc20,lw20,iwrkd,ldiwrk,ier) c nur=0 nuz=2 call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv, . nz,ffvpsi,cc02,lw02,iwrkd,ldiwrk,ier) c nur=1 nuz=1 call dierckx_parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv, . nz,ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier) c c scaling of f_poloidal c do i=1,nrho fpol(i)=fpol(i)*factb wf(i)=1.0d0 end do wf(nrho)=1.0d2 c c spline interpolation of fpol(i) c iopt=0 xb=0.0d0 xe=1.0d0 ssfp=0.01d0 call dierckx_curfit(iopt,nrho,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest, . nsft,tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier) c call dierckx_splev(tfp,nsft,cfp,3,psinr,fpoli,nrho,ier) fpolas=fpoli(nrho) c c no limiter shape provided yet nlim=0 c c compute max and min z of last closed surface c rbmin=rmaxis rbmax=rmaxis if (nbnd.gt.1) then zbmin=1.0d+30 zbmax=-1.0d+30 do i=1,nbnd if(zbnd(i).le.zbmin) then zbmin=zbnd(i) rbmin=rbnd(i) end if if(zbnd(i).ge.zbmax) then zbmax=zbnd(i) rbmax=rbnd(i) end if end do else zbmin=-1.0d+30 zbmax=1.0d+30 end if if(zbmin.le.zmnm) zbmin=zbmin+dz if(rbmin.le.rmnm) rbmin=rbmin+dr if(zbmax.ge.zmxm) zbmax=zbmax-dz if(rbmax.ge.rmxm) rbmax=rbmax-dr c c start with uncorrected normalized psi c psinop=0.0d0 psinxp=1.0d0 psiant=1.0d0 c c search for O-point c call points_ox(rmaxis,zmaxis,rmop,zmop,psinoptmp,info) rmaxis=rmop zmaxis=zmop write(*,'(a,2f8.4,es12.5)') 'O-point',rmop,zmop,psinoptmp c c search for X-point if ixp not = 0 c if(ixp.ne.0) then if(ixp.lt.0) then r10=rbmin z10=zbmin call points_ox(r10,z10,rxp,zxp,psinxptmp,info) if(psinxp.ne.-1.0d0) then write(*,'(a,2f8.4,es12.5)') 'X-point',rxp,zxp,psinxptmp rbmin=rxp zbmin=zxp psinop=psinoptmp psinxp=psinxptmp psiant=psinxp-psinop psin1=1.0d0 r10=rmaxis z10=(zbmax+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) rbmax=r1 zbmax=z1 else ixp=0 c write(*,'(a)') 'no X-point' end if else r10=rmop z10=zbmax call points_ox(r10,z10,rxp,zxp,psinxptmp,info) if(psinxp.ne.-1.0d0) then write(*,'(a,2f8.4,e16.8)') 'X-point',rxp,zxp,psinxptmp zbmax=zxp psinop=psinoptmp psinxp=psinxptmp psiant=psinxp-psinop psin1=1.0d0 z10=(zbmin+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) zbmin=z1 else ixp=0 c write(*,'(a)') 'no X-point' end if end if end if c if (ixp.eq.0) then psin1=1.0d0 psinop=psinoptmp psiant=psin1-psinop r10=rmaxis z10=(zbmax+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) zbmax=z1 rbmax=r1 z10=(zbmin+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) zbmin=z1 rbmin=r1 write(*,'(a,4f8.4)') 'no X-point ',rbmin,zbmin,rbmax,zbmax end if write(*,*) ' ' c compute B_toroidal on axis btaxis=fpol(1)/rmaxis btrcen=btrcen*factb write(*,'(a,f8.4)') 'factb = ',factb write(*,'(a,f8.4)') 'BT_centr= ',btrcen write(*,'(a,f8.4)') 'BT_axis = ',btaxis c compute normalized rho_tor from eqdsk q profile call rhotor(nrho) phitedge=abs(psia)*rhotsx*2*pi rrtor=sqrt(abs(phitedge/btrcen)/pi) call rhopol c write(*,*) rhotsx,phitedge,rrtor,abs(psia) c compute flux surface averaged quantities rup=rmaxis rlw=rmaxis zup=zmaxis+(zbmax-zmaxis)/10.0d0 zlw=zmaxis-(zmaxis-zbmin)/10.0d0 call flux_average c ipr=1 c call contours_psi(1.0d0,np,rcon,zcon,ipr) c do ii=1,2*np+1 c write(632,*) rcon(ii), zcon(ii) c end do c c locate psi surface for q=1.5 and q=2 rup=rmaxis rlw=rmaxis zup=(zbmax+zmaxis)/2.0d0 zlw=(zmaxis+zbmin)/2.0d0 q2=2.0d0 q15=1.5d0 call vmaxmini(qpsi,nrho,qmin,qmax,iqmn,iqmx) if (q15.gt.qmin.and.q15.lt.qmax) then call surfq(q15,psi15) rhot15=frhotor(psi15) write(*,'(3(a,f8.5))') 'psi_1.5 = ',psi15,' rhop_1.5 = ' . ,sqrt(psi15),' rhot_1.5 = ',rhot15 end if if (q2.gt.qmin.and.q2.lt.qmax) then call surfq(q2,psi2) rhot2=frhotor(psi2) write(*,'(3(a,f8.5))') 'psi_2 = ',psi2,' rhop_2 = ' . ,sqrt(psi2),' rhot_2 = ',rhot2 end if c c locate btot=bres c call bfield_res c return end c c c subroutine points_ox(rz,zz,rf,zf,psinvf,info) use const_and_precisions, only : comp_eps implicit real*8 (a-h,o-z) parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) external fcnox common/psinv/psinv xvec(1)=rz xvec(2)=zz tol = sqrt(comp_eps) call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) if(info.gt.1) then write(*,'(a,i2,a,2f8.4)') ' info subr points_ox =',info, . ' O/X coord.',xvec end if rf=xvec(1) zf=xvec(2) psinvf=psinv return end c c c subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/pareq1/psia call equinum(x(1),x(2)) if (iflag.ne.2) then fvec(1) = dpsidr/psia fvec(2) = dpsidz/psia else fjac(1,1) = ddpsidrr/psia fjac(1,2) = ddpsidrz/psia fjac(2,1) = ddpsidrz/psia fjac(2,2) = ddpsidzz/psia end if return end c c c subroutine points_tgo(rz,zz,rf,zf,psin,info) use const_and_precisions, only : comp_eps implicit real*8 (a-h,o-z) parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) external fcntgo common/cnpsi/h h=psin xvec(1)=rz xvec(2)=zz tol = sqrt(comp_eps) call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) if(info.gt.1) then end if rf=xvec(1) zf=xvec(2) return end c c c subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag) implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/psinv/psinv common/cnpsi/h common/pareq1/psia call equinum(x(1),x(2)) if (iflag.ne.2) then fvec(1) = psinv-h fvec(2) = dpsidr/psia else fjac(1,1) = dpsidr/psia fjac(1,2) = dpsidz/psia fjac(2,1) = ddpsidrr/psia fjac(2,2) = ddpsidrz/psia end if return end c c c subroutine print_prof implicit real*8 (a-h,o-z) parameter(nnw=501,eps=1.d-4) dimension psinr(nnw),qpsi(nnw) c common/psinr/psinr common/qpsi/qpsi common/eqnn/nr,nz,nrho,npp,nintp common/dddens/dens,ddens c write(645,*) ' #psi rhot ne Te q Jphi' psin=0.0d0 rhop=0.0d0 rhot=0.0d0 call density(psin) call tor_curr_psi(eps,ajphi) te=temperature(psin) qq=qpsi(1) c write(645,111) psin,rhot,dens,te,qq,ajphi*1.d-6 c nst=nrho do i=2,nst psin=dble(i-1)/dble(nst-1) rhop=sqrt(psin) c call density(psin) te=temperature(psin) c call vlocate(psinr,nrho,psin,ips) ips=min(max(ips,1),nrho-1) call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1), . psin,qq) rhot=frhotor(psin) call tor_curr_psi(psin,ajphi) write(645,111) psin,rhot,dens,te,qq,ajphi*1.d-6 end do c return 111 format(12(1x,e12.5)) end c c c subroutine surfq(qval,psival) implicit real*8 (a-h,o-z) parameter(nnw=501) parameter(ncnt=100,ncntt=2*ncnt+1) dimension psinr(nnw),qpsi(nnw) dimension rcon(ncntt),zcon(ncntt) c common/psinr/psinr common/qpsi/qpsi common/eqnn/nr,nz,nrho,npp,nintp c c locate psi surface for q=qval c call vlocate(qpsi,nrho,qval,i1) call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1), . qval,psival) ipr=1 call contours_psi(psival,ncnt,rcon,zcon,ipr) return end c c c subroutine bfield_res implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) dimension rv(nnw),zv(nnh),psin(nnw,nnh) dimension btotal(nnw,nnh) parameter(icmx=2002) dimension rrcb(icmx),zzcb(icmx),ncpts(10) c common/cpsin/rv,zv,psin common/eqnn/nr,nz,nrho,npp,nintp common/parbres/bres common/btt/btotal c c Btotal on psi grid c btmx=-1.0d30 btmn=1.0d30 do j=1,nr rrj=rv(j) do k=1,nz zzk=zv(k) call bfield(rrj,zzk,bbphi,bbr,bbz) btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2) if(btotal(j,k).ge.btmx) btmx=btotal(j,k) if(btotal(j,k).le.btmn) btmn=btotal(j,k) c write(640,113) j,rrj,zzk,btotal(j,k) enddo c write(640,*) ' ' enddo c c compute Btot=Bres and Btot=Bres/2 c write(630,*)'#i Btot R z' do n=1,5 bbb=bres/dble(n) if (bbb.ge.btmn.and.bbb.le.btmx) then call cniteq(bbb,nconts,ncpts,nctot,rrcb,zzcb,1) do inc=1,nctot write(630,113) inc,bbb,rrcb(inc),zzcb(inc) end do end if write(630,*) ' ' end do c return 113 format(i6,12(1x,e12.5)) end c c subroutine profdata(mrho, psijet, te, dne, zeff) implicit real*8 (a-h,o-z) c integer mrho real*8 psijet(mrho), te(mrho), dne(mrho), zeff(mrho) c parameter(npmx=501,npest=npmx+4) dimension psrad(npmx),terad(npmx),derad(npmx),zfc(npmx) dimension ct(npmx,4),cz(npmx,4) parameter(lwrkf=npmx*4+npest*16) dimension tfn(npest),cfn(npest),wrkf(lwrkf),iwrkf(npest),wf(npmx) dimension densi(npest),ddensi(npest),d2densi(npest),wrkfd(npest) c common/densbnd/psdbnd common/denspp/psnpp,aa,bb,cc common/eqnn/nr,nz,nrho,npp,nintp common/iipr/iprof common/nfile/neqdsk,nprof common/crad/psrad,derad,terad,zfc common/coffte/ct common/coffz/cz common/factb/factb common/coffdt/tfn common/coffdnst/nsfd common/cofffn/cfn common/scal/iscal common/facttn/factt,factn c c read plasma profiles from file if iprof>0 c if (iscal.eq.0) then aat=2.0d0/3.0d0 aan=4.0d0/3.0d0 else aat=1.0d0 aan=1.0d0 end if ffact=factb if(iscal.eq.2) ffact=1.0d0 c if (iprof.gt.0) then do i=1,mrho psrad(i)=psijet(i) terad(i)=te(i)*1.d-3*ffact**aat*factt derad(i)=dne(i)*1.d-19*ffact**aan*factn zfc(i)=zeff(i) wf(i)=1.0d0 end do psrad(1)=0.0d0 npp=mrho c c spline approximation of temperature and Zeff c iopt=0 call difcsn(psrad,terad,npmx,npp,iopt,ct,ier) c iopt=0 call difcsn(psrad,zfc,npmx,npp,iopt,cz,ier) c c spline approximation of density c iopt=0 xb=0.0d0 xe=psrad(npp) kspl=3 sspl=.001d0 c call dierckx_curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest, . nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier) c call dierckx_splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier) nu=1 call dierckx_splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier) dnpp=densi(npp) ddnpp=ddensi(npp) c nu=2 call dierckx_splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier) d2dnpp=d2densi(npp) if(derad(npp).eq.0.0d0) then psdbnd=psrad(npp) else psnpp=psrad(npp) dpsb=-psnpp+psdbnd nn=3 nn1=nn+1 nn2=nn+2 aa=(nn1*nn2*dnpp+2*nn1*ddnpp*dpsb+d2dnpp*dpsb**2) aa=aa/(-dpsb)**nn/2.0d0 bb=-(nn*nn2*dnpp+(2*nn+1)*ddnpp*dpsb+d2dnpp*dpsb**2) bb=bb/(-dpsb)**nn1 cc=(nn1*nn*dnpp+2*nn*ddnpp*dpsb+d2dnpp*dpsb**2) cc=cc/(-dpsb)**nn2/2.0d0 end if c end if c return end c c c subroutine rhotor(nr) implicit real*8(a-h,o-z) parameter(nnw=501) dimension psinr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),crhot(nnw,4) common/psinr/psinr common/qpsi/qpsi common/rhotsx/rhotsx common/crhot/crhot common/cq/cq c c normalized toroidal rho : ~ Integral q(psi) dpsi c iopt=0 call difcsn(psinr,qpsi,nnw,nr,iopt,cq,ier) c rhotnr(1)=0.0d0 do k=1,nr-1 dx=psinr(k+1)-psinr(k) drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0d0+dx*(cq(k,3)/3.0d0 . +dx*cq(k,4)/4.0d0))) rhotnr(k+1)=rhotnr(k)+drhot end do rhotsx=rhotnr(nr) do k=2,nr rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr)) end do c c spline interpolation of rhotor c iopt=0 call difcsn(psinr,rhotnr,nnw,nr,iopt,crhot,ier) return end function fq_eq(psi) implicit real*8(a-h,o-z) parameter(nnw=501) dimension psinr(nnw),cq(nnw,4) common/psinr/psinr common/eqnn/nr,nz,nrho,npp,nintp common/cq/cq call vlocate(psinr,nrho,psi,irt) irt=max(1,min(irt,nrho-1)) dps=psi-psinr(irt) fq_eq=fspli(cq,nnw,irt,dps) return end function frhotor_eq(psi) implicit real*8(a-h,o-z) parameter(nnw=501) dimension psinr(nnw),crhot(nnw,4) common/psinr/psinr common/eqnn/nr,nz,nrho,npp,nintp common/crhot/crhot c call vlocate(psinr,nrho,psi,irt) irt=max(1,min(irt,nrho-1)) dps=psi-psinr(irt) frhotor_eq=fspli(crhot,nnw,irt,dps) return end function frhotor(psi) implicit real*8(a-h,o-z) frhotor=frhotor_eq(psi) return end function frhotor_av(psi) implicit real*8(a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),crhotq(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,nrho,npp,nintp common/crhotq/crhotq rpsi=sqrt(psi) c ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 call vlocate(rpstab,nintp,rpsi,ip) ip=max(1,min(ip,nintp-1)) dps=rpsi-rpstab(ip) frhotor_av=fspli(crhotq,nnintp,ip,dps) return end subroutine rhopol implicit real*8(a-h,o-z) parameter(nnr=101,nrest=nnr+4) parameter(lwrkp=nnr*4+nrest*16) dimension rhop(nnr),rhot(nnr),rhopi(nnr) dimension trp(nrest),crp(nrest),wp(nrest) dimension wrkp(lwrkp),iwrkp(nrest) common/coffrtp/trp common/coffrn/nsrp common/coffrp/crp dr=1.0d0/dble(nnr-1) do i=2,nnr-1 rhop(i)=(i-1)*dr psin=rhop(i)*rhop(i) rhot(i)=frhotor_eq(psin) wp(i)=1.0d0 end do rhop(1)=0.0d0 rhot(1)=0.0d0 wp(1)=1.0d3 rhop(nnr)=1.0d0 rhot(nnr)=1.0d0 wp(nnr)=1.0d3 c spline interpolation of rhopol versus rhotor iopt=0 xb=0.0d0 xe=1.0d0 ss=0.00001 kspl=3 call dierckx_curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest, . nsrp,trp,crp,rp,wrkp,lwrkp,iwrkp,ier) c write(*,*) ier call dierckx_splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier) do i=1,nnr write(644,*) rhop(i),rhot(i),rhopi(i) end do return end function frhopol(rhot) implicit real*8(a-h,o-z) parameter(nnr=101,nrest=nnr+4) dimension trp(nrest),crp(nrest),rrs(1),ffspl(1) common/coffrtp/trp common/coffrn/nsrp common/coffrp/crp rrs(1)=rhot call dierckx_splev(trp,nsrp,crp,3,rrs,ffspl,1,ier) frhopol=ffspl(1) return end subroutine cniteq(h, ncon, npts, icount, rcon, zcon,ichoi) implicit real*8 (a-h,o-z) c c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. c (based on an older code) c parameter(nnw=501,nnh=501,icmx=2002,nna=nnw*nnh) dimension a(nna),ja(3,2),lx(1000),npts(10) dimension rcon(icmx),zcon(icmx) dimension rqgrid(nnw),zqgrid(nnh),psin(nnw,nnh),btotal(nnw,nnh) c common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/cpsin/rqgrid,zqgrid,psin common/btt/btotal common/eqnn/nr,nz,nrho,npp,nintp c data px/0.5d0/ c if(ichoi.eq.0) then do j=1,nz do i=1,nr a(nr*(j-1)+i)=psin(i,j) enddo enddo endif c if(ichoi.eq.1) then do j=1,nz do i=1,nr a(nr*(j-1)+i)=btotal(i,j) enddo enddo endif c do ico=1,icmx rcon(ico)=0.0d0 zcon(ico)=0.0d0 enddo c nrqmax=nr nreq=nr nzeq=nz drgrd=dr dzgrd=dz c ncon = 0 do i=1,10 npts(i) = 0 end do iclast = 0 icount = 0 mpl=0 ix=0 mxr = nrqmax * (nzeq - 1) n1 = nreq - 1 c do 3 jx=2,n1 do 3 jm=jx,mxr,nrqmax j = jm + nrqmax ah=a(j)-h if(ah) 60,60,61 60 if(a(jm)-h) 62,62,1 61 if(a(jm)-h) 1,1,63 1 ix=ix+1 lx(ix)=-j if(ah) 62,62,63 62 if(a(j-1)-h) 3,3,4 63 if(a(j-1)-h) 4,4,3 4 ix=ix+1 lx(ix)=j 3 continue c do 79 jm=nreq,mxr,nrqmax j = jm + nrqmax ah=a(j)-h if(ah) 64,64,65 64 if(a(j-1)-h) 66,66,76 65 if(a(j-1)-h) 76,76,67 76 ix=ix+1 lx(ix)=j if(ah) 66,66,67 66 if(a(jm)-h) 79,79,78 67 if(a(jm)-h) 78,78,79 78 ix=ix+1 lx(ix)=-j 79 continue c do 75 jm=1,mxr,nrqmax j = jm + nrqmax if(a(j)-h) 68,68,69 68 if(a(jm)-h)75,75,71 69 if(a(jm)-h)71,71,75 71 ix=ix+1 lx(ix) =-j 75 continue c do 73 j=2,nreq if(a(j)-h) 58,58,59 58 if(a(j-1)-h) 73,73,72 59 if(a(j-1)-h) 72,72,73 72 ix=ix+1 lx(ix)=j 73 continue c if(ix) 50,50,8 108 if(mpl.lt.4) go to 8 ncon = ncon + 1 npts(ncon) = icount - iclast iclast = icount 8 in=ix jx=lx(in) jfor=0 lda=1 ldb=2 30 if(jx) 21,22,22 21 jabs=-jx jnb = jabs - nrqmax go to 23 22 jabs=jx jnb=jabs-1 23 adn=a(jabs)-a(jnb) if(adn) 24,9,24 24 px=(a(jabs)-h)/adn 9 kx = (jabs - 1) / nrqmax ikx = jabs - nrqmax * kx - 1 if(jx) 25,26,26 25 x = drgrd * ikx y = dzgrd * (kx - px) go to 27 26 x = drgrd * (ikx - px) y = dzgrd * kx 27 continue icount = icount + 1 rcon(icount) = x + rqgrid(1) zcon(icount) = y + zqgrid(1) mpl= icount itm = 1 ja(1,1) = jabs + nrqmax j=1 if(jx) 10,10,11 10 ja(1,1) = -jabs-1 j=2 11 ja(2,1)=-ja(1,1) ja(3,1) = -jx + 1 - nrqmax ja(3,2) = -jx ja(j,2) = jabs - nrqmax k= 3-j ja(k,2) = 1-jabs if (kx) 14,14,39 39 if(ikx) 14,14,36 36 if(ikx + 1 - nreq) 35, 37, 37 37 if(jx) 38,38,35 35 if(jfor) 28,29,28 28 do 13 i=1,3 if(jfor-ja(i,2)) 13,14,13 13 continue 38 lda=2 go to 15 14 lda=1 15 ldb=lda 29 do 18 k=1,3 do 18 l=lda,ldb do 16 i=1,ix if(lx(i)-ja(k,l)) 16,17,16 16 continue go to 18 17 itm=itm+1 inext= i if(jfor) 19,33,19 33 if(itm .gt. 3) goto 20 18 continue 19 lx(in)=0 if(itm .eq. 1) goto 6 20 jfor=jx jx=lx(inext) in = inext go to 30 6 if(lx(ix)) 108,7,108 7 ix= ix-1 if(ix) 51,51,6 51 if(mpl.lt.4) go to 50 ncon = ncon + 1 npts(ncon) = icount - iclast iclast = icount 50 continue c return end c c c subroutine contours_psi(h,np,rcn,zcn,ipr) use const_and_precisions, only : pi implicit real*8 (a-h,o-z) parameter(mest=4,kspl=3) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) dimension rcn(2*np+1),zcn(2*np+1) dimension cc(nnw*nnh),tr(nrest),tz(nzest) dimension czc(nrest),zeroc(mest) c common/pareq1/psia common/pareq1t/psiant,psinop common/coffeqn/nsr,nsz,nsft common/coffeq/cc common/coffeqt/tr,tz common/cnt/rup,zup,rlw,zlw common/rwallm/rwallm c ra=rup rb=rlw za=zup zb=zlw call points_tgo(ra,za,rup,zup,h,info) call points_tgo(rb,zb,rlw,zlw,h,info) th=pi/dble(np) rcn(1)=rlw zcn(1)=zlw rcn(2*np+1)=rlw zcn(2*np+1)=zlw rcn(np+1)=rup zcn(np+1)=zup do ic=2,np zc=zlw+(zup-zlw)*(1.0d0-cos(th*(ic-1)))/2.0d0 iopt=1 call dierckx_sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc, . ier) if(ier.gt.0) write(*,*) ' dierckx_sprofil =',ier val=h*psiant+psinop call dierckx_sproota(val,tr,nsr,czc,zeroc,mest,m,ier) if (zeroc(1).gt.rwallm) then rcn(ic)=zeroc(1) zcn(ic)=zc rcn(2*np+2-ic)=zeroc(2) zcn(2*np+2-ic)=zc else rcn(ic)=zeroc(2) zcn(ic)=zc rcn(2*np+2-ic)=zeroc(3) zcn(2*np+2-ic)=zc end if end do if (ipr.gt.0) then do ii=1,2*np+1 write(631,111) ii,h,rcn(ii),zcn(ii) end do write(631,*) write(631,*) end if return 111 format(i6,12(1x,e12.5)) 99 format(12(1x,e12.5)) end c c c subroutine flux_average use const_and_precisions, only : pi, mu0=>mu0_ implicit real*8 (a-h,o-z) real*8 lam parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41) parameter(zero=0.0d0,one=1.0d0) parameter(ccj=1.0d0/mu0) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ . njest*nnintp+nlest+54) parameter(kwrk=nnintp+nlam+njest+nlest+3) parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp) dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp) dimension dadrhotv(nnintp),cdadrhot(nnintp,4) dimension dvdrhotv(nnintp),cdvdrhot(nnintp,4) dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp) dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp) dimension rcon(2*ncnt+1),zcon(2*ncnt+1) dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1) dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4) dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4) dimension cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4) dimension cfc(nnintp,4),crhotq(nnintp,4) dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp) dimension alam(nlam),fhlam(nnintp,nlam) dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam) dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)) dimension iwrk(kwrk),wrk(lwrk) dimension ch01(lw01),weights(nlam) c common/cent/btrcen,rcen common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/eqnn/nr,nz,nrho,npp,nintp common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv c common/pstab/rpstab common/flav/vvol,rri,rbav,bmxpsi,bmnpsi common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc common/cratj/cratja,cratjb,cratjpl common/crhotq/crhotq common/cnt/rup,zup,rlw,zlw common/zbound/zbmin,zbmax common/rarea/rarea common/phitedge/phitedge common/cdadrhot/cdadrhot common/cdvdrhot/cdvdrhot c common/coffvl/ch common/coffdvl/ch01 common/coffvlt/tjp,tlm common/coffvln/njpt,nlmt c computation of flux surface averaged quantities write(631,*)' #i psin R z' nintp=nnintp ninpr=(nintp-1)/10 dlam=1.0d0/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam fhlam(1,l)=sqrt(1.0d0-alam(l)) ffhlam(l)=fhlam(1,l) dffhlam(l)=-0.5d0/sqrt(1.0d0-alam(l)) weights(l)=1.0d0 end do weights(1)=0.5d0 weights(nlam)=0.5d0 alam(nlam)=1.0d0 fhlam(1,nlam)=0.0d0 ffhlam(nlam)=0.0d0 dffhlam(nlam)=-99999.0d0 jp=1 anorm=2.0d0*pi*rmaxis/abs(btaxis) b2av=btaxis**2 dvdpsi=2.0d0*pi*anorm dadpsi=2.0d0*pi/abs(btaxis) ratio_cdator=abs(btaxis/btrcen) ratio_cdbtor=1.0d0 ratio_pltor=1.0d0 qq=1.0d0 fc=1.0d0 pstab(1)=0.0d0 rpstab(1)=0.0d0 rhotqv(1)=0.0d0 vcurrp(1)=0.0d0 vajphiav(1)=0.0d0 bmxpsi(1)=abs(btaxis) bmnpsi(1)=abs(btaxis) bav(1)=abs(btaxis) rbav(1)=1.0d0 rri(1)=rmaxis varea(1)=0.0d0 vvol(1)=0.0d0 vratjpl(1)=ratio_pltor vratja(1)=ratio_cdator vratjb(1)=ratio_cdbtor ffc(1)=fc rhot2q=0.0d0 C rup=rmaxis C rlw=rmaxis C zup=(zbmax+zmaxis)/2.0d0 C zlw=(zmaxis+zbmin)/2.0d0 do jp=2,nintp height=dble(jp-1)/dble(nintp-1) if(jp.eq.nintp) height=0.9999d0 ipr=0 jpr=mod(jp,ninpr) if(jpr.eq.1) ipr=1 height2=height*height call contours_psi(height2,ncnt,rcon,zcon,ipr) c r2iav=0.0d0 anorm=0.0d0 dadpsi=0.0d0 currp=0.0d0 b2av=0.0d0 area=0.0d0 volume=0.0d0 ajphiav=0.0d0 bbav=0.0d0 bmmx=-1.0d+30 bmmn=1.0d+30 call bfield(rcon(1),zcon(1),bphi,brr,bzz) call tor_curr(rcon(1),zcon(1),ajphi0) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) bv(1)=btot0 bpv(1)=bpoloid0 rpsim0=rcon(1) do inc=1,ncntt-1 inc1=inc+1 dla=sqrt((rcon(inc)-rmaxis)**2+(zcon(inc)-zmaxis)**2) dlb=sqrt((rcon(inc1)-rmaxis)**2+(zcon(inc1)-zmaxis)**2) dlp=sqrt((rcon(inc1)-rcon(inc))**2+ . (zcon(inc1)-zcon(inc))**2) drc=(rcon(inc1)-rcon(inc)) c compute length, area and volume defined by psi=height^2 ph=0.5d0*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) rzp=rcon(inc1)*zcon(inc1) rz=rcon(inc)*zcon(inc) volume=pi*(rzp+rz)*drc+volume c compute line integral on the contour psi=height^2 rpsim=rcon(inc1) zpsim=zcon(inc1) call bfield(rpsim,zpsim,bphi,brr,bzz) call tor_curr(rpsim,zpsim,ajphi) btot=sqrt(bphi**2+brr**2+bzz**2) bpoloid=sqrt(brr**2+bzz**2) dlpv(inc)=dlp bv(inc1)=btot bpv(inc1)=bpoloid dlph=0.5d0*dlp anorm=anorm+dlph*(1.0d0/bpoloid+1.0d0/bpoloid0) dadpsi=dadpsi+dlph* . (1.0d0/(bpoloid*rpsim)+1.0d0/(bpoloid0*rpsim0)) currp=currp+dlph*(bpoloid+bpoloid0) b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) r2iav=r2iav+dlph* . (1.0d0/(bpoloid*rpsim**2)+1.0d0/(bpoloid0*rpsim0**2)) ajphiav=ajphiav+dlph* . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) ajphi0=ajphi rpsim0=rpsim bpoloid0=bpoloid btot0=btot c computation maximum/minimum B values on given flux surface if(btot.le.bmmn) bmmn=btot if(btot.ge.bmmx) bmmx=btot end do c bav= [T] , b2av= [T^2] , rbav=/b_min c anorm = int d l_p/B_p = dV/dpsi/(2pi) c r2iav=<1/R^2> [m^-2] , c riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), c rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] c currp = plasma current within psi=const bbav=bbav/anorm r2iav=r2iav/anorm dvdpsi=2.0d0*pi*anorm riav=dadpsi/anorm b2av=b2av/anorm vcurrp(jp)=ccj*currp vajphiav(jp)=ajphiav/dadpsi c area == varea, volume == vvol c flux surface minor radius == (area/pi)^1/2 c ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 c ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / c ratio_pltor = Jcd_||/J_phi Jcd_|| = pstab(jp)=height2 rpstab(jp)=height vvol(jp)=abs(volume) varea(jp)=area bav(jp)=bbav rbav(jp)=bbav/bmmn bmxpsi(jp)=bmmx bmnpsi(jp)=bmmn rri(jp)=bav(jp)/abs(fpolv*r2iav) ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) vratjpl(jp)=ratio_pltor vratja(jp)=ratio_cdator vratjb(jp)=ratio_cdbtor qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) qqv(jp)=qq dadrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) . *dadpsi/pi dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) . *dvdpsi/pi c c write(647,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp) c computation of rhot from calculated q profile rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) . /dble(nintp-1) rhotqv(jp)=sqrt(rhot2q) c computation of fraction of circulating/trapped fraction fc, ft c and of function H(lambda,rhop) c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ fc=0.0d0 shlam=0.0d0 do l=nlam,1,-1 lam=alam(l) srl=0.0d0 rl2=1.0d0-lam*bv(1)/bmmx rl0=0.d0 if(rl2.gt.0) rl0=sqrt(rl2) do inc=1,ncntt-1 rl2=1.0d0-lam*bv(inc+1)/bmmx rl=0.0d0 if(rl2.gt.0) rl=sqrt(rl2) srl=srl+0.5d0*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) rl0=rl end do srl=srl/anorm dhlam=0.5d0/srl fc=fc+lam/srl*weights(l) if(l.eq.nlam) then fhlam(jp,l)=0.0d0 ffhlam(nlam*(jp-1)+l)=0.0d0 dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam else shlam=shlam+0.5d0*(dhlam+dhlam0)*dlam fhlam(jp,l)=shlam dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam end if end do fc=0.75d0*b2av/bmmx**2*fc*dlam ffc(jp)=fc ccfh=bmmn/bmmx/fc do l=1,nlam ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) end do end do write(646,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// .' Area Vol |I_pl| qq fc ratioJa ratioJb' qqv(1)=qqv(2) vajphiav(1)=vajphiav(2) do jp=1,nintp rhotqv(jp)=rhotqv(jp)/rhotqv(nintp) if(jp.eq.nintp) then rhotqv(jp)=1.0d0 rpstab(jp)=1.0d0 pstab(jp)=1.0d0 end if rhot_eq=frhotor_eq(pstab(jp)) write(646,99) pstab(jp),rhot_eq,rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp) . ,vratja(jp),vratjb(jp) end do rarea=sqrt(varea(nintp)/pi) c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi c used for computations of dP/dV and J_cd c spline coefficients of rhot iopt=0 call difcsg(rpstab,vvol,nintp,iopt,cvol,ier) iopt=0 call difcsg(rpstab,rbav,nintp,iopt,crbav,ier) iopt=0 call difcsg(rpstab,rri,nintp,iopt,crri,ier) iopt=0 call difcsg(rpstab,bmxpsi,nintp,iopt,cbmx,ier) iopt=0 call difcsg(rpstab,bmnpsi,nintp,iopt,cbmn,ier) iopt=0 call difcsg(rpstab,vratja,nintp,iopt,cratja,ier) iopt=0 call difcsg(rpstab,vratjb,nintp,iopt,cratjb,ier) iopt=0 call difcsg(rpstab,vratjpl,nintp,iopt,cratjpl,ier) iopt=0 call difcsg(rpstab,varea,nintp,iopt,carea,ier) iopt=0 call difcsg(rpstab,ffc,nintp,iopt,cfc,ier) iopt=0 call difcsg(rpstab,rhotqv,nintp,iopt,crhotq,ier) iopt=0 call difcsg(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier) iopt=0 call difcsg(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,ier) c spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 s=0.0d0 call dierckx_regrid(iopt,nintp,rpstab,nlam,alam,ffhlam, . zero,one,zero,one,ksp,ksp,s, . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) njpt=njp nlmt=nlm call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1, . ch01,lw01,ier) return 99 format(20(1x,e12.5)) end function fdadrhot(rpsi) implicit real*8(a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),cdadrhot(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,nrho,npp,nintp common/cdadrhot/cdadrhot c ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 call vlocate(rpstab,nintp,rpsi,ip) ip=max(1,min(ip,nintp-1)) dps=rpsi-rpstab(ip) fdadrhot=fspli(cdadrhot,nnintp,ip,dps) return end function fdvdrhot(rpsi) implicit real*8(a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),cdvdrhot(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,nrho,npp,nintp common/cdvdrhot/cdvdrhot c ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 call vlocate(rpstab,nintp,rpsi,ip) ip=max(1,min(ip,nintp-1)) dps=rpsi-rpstab(ip) fdvdrhot=fspli(cdvdrhot,nnintp,ip,dps) return end subroutine vectinit use beamdata, only : jmx,kmx,nmx,psjki,tauv,alphav,pdjki,ppabs, . currj,didst,ccci,iiv,iop,iow,tau1v, . nrayr,nrayth,nstep,alloc_beam implicit real*8 (a-h,o-z) parameter(tmax=5,npts=500) real*8 ttv(npts+1),extv(npts+1) common/ttex/ttv,extv c common/warm/iwarm,ilarm c common/ierr/ierr c if(nstep.gt.nmx) nstep=nmx if(nrayr.gt.jmx) nrayr=jmx if(nrayth.gt.kmx) nrayth=kmx call alloc_beam(ierr) if (ierr.ne.0) return c do i=1,nstep do k=1,nrayth do j=1,nrayr psjki(j,k,i)=0.0d0 tauv(j,k,i)=0.0d0 alphav(j,k,i)=0.0d0 pdjki(j,k,i)=0.0d0 ppabs(j,k,i)=0.0d0 didst(j,k,i)=0.0d0 ccci(j,k,i)=0.0d0 currj(j,k,i)=0.0d0 iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 tau1v(j,k)=0.0d0 end do end do end do c if(iwarm.gt.1) then dt=2.0d0*tmax/dble(npts) do i = 1, npts+1 ttv(i) = -tmax+dble(i-1)*dt extv(i)=exp(-ttv(i)*ttv(i)) end do end if c return end subroutine vectfree use beamdata, only : dealloc_beam implicit none c call dealloc_beam c return end subroutine vectinit2 use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj, . didst,ccci,iiv,iop,iow,nrayr,nrayth,nstep implicit real*8 (a-h,o-z) do i=1,nstep do k=1,nrayth do j=1,nrayr psjki(j,k,i)=0.0d0 tauv(j,k,i)=0.0d0 alphav(j,k,i)=0.0d0 pdjki(j,k,i)=0.0d0 ppabs(j,k,i)=0.0d0 didst(j,k,i)=0.0d0 ccci(j,k,i)=0.0d0 currj(j,k,i)=0.0d0 iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 end do end do end do return end c c c subroutine paraminit implicit real*8(a-h,o-z) c common/istep/istep common/istgr/istpr,istpl common/ierr/ierr common/istop/istop common/ipol/ipolc c istpr=0 istpl=1 ierr=0 istep=0 istop=0 ipolc=0 c return end c c subroutine updatepos use beamdata, only : xc,xco,du1,du1o,ywrk,nrayr,nrayth c implicit real*8 (a-h,o-z) c do j=1,nrayr do k=1,nrayth if(j.eq.1.and.k.gt.1) then xco(1,j,k)=xco(1,j,1) xco(2,j,k)=xco(2,j,1) xco(3,j,k)=xco(3,j,1) xc(1,j,k)=xc(1,j,1) xc(2,j,k)=xc(2,j,1) xc(3,j,k)=xc(3,j,1) else xco(1,j,k)=xc(1,j,k) xco(2,j,k)=xc(2,j,k) xco(3,j,k)=xc(3,j,k) xc(1,j,k)=ywrk(1,j,k) xc(2,j,k)=ywrk(2,j,k) xc(3,j,k)=ywrk(3,j,k) end if du1o(1,j,k)=du1(1,j,k) du1o(2,j,k)=du1(2,j,k) du1o(3,j,k)=du1(3,j,k) end do end do c return end c c c subroutine gradi use beamdata, only : dffiu,ddffiu,xc,xco,du1,du1o,gri,ggri, . dgrad2v,grad2,nrayr,nrayth implicit real*8 (a-h,o-z) dimension dxv1(3),dxv2(3),dxv3(3),dgu(3) dimension dgg1(3),dgg2(3),dgg3(3) dimension df1(3),df2(3),df3(3) c c compute grad u and grad(S_I) c do k=1,nrayth do j=1,nrayr if(j.eq.1) then gri(1,j,k)=0.0d0 gri(2,j,k)=0.0d0 gri(3,j,k)=0.0d0 jp=j+1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,jp,k)-xc(iv,j,k) dxv2(iv)=xc(iv,jp,kp)-xc(iv,jp,km) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do call solg0(dxv1,dxv2,dxv3,dgu) else jm=j-1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do call solg0(dxv1,dxv2,dxv3,dgu) end if du1(1,j,k)=dgu(1) du1(2,j,k)=dgu(2) du1(3,j,k)=dgu(3) gri(1,j,k)=dgu(1)*dffiu(j) gri(2,j,k)=dgu(2)*dffiu(j) gri(3,j,k)=dgu(3)*dffiu(j) grad2(j,k)=gri(1,j,k)**2+gri(2,j,k)**2+gri(3,j,k)**2 end do end do c c compute derivatives of grad u and grad(S_I) c do k=1,nrayth do j=1,nrayr if(j.eq.1) then jp=j+1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,jp,kp)-xc(iv,jp,km) dxv2(iv)=xc(iv,jp,k)-xc(iv,j,k) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do df1(1)=du1(1,jp,kp)-du1(1,jp,km) df1(2)=du1(1,jp,k)-du1(1,j,k) df1(3)=du1(1,j,k)-du1o(1,j,k) df2(1)=du1(2,jp,kp)-du1(2,jp,km) df2(2)=du1(2,jp,k)-du1(2,j,k) df2(3)=du1(2,j,k)-du1o(2,j,k) df3(1)=du1(3,jp,kp)-du1(3,jp,km) df3(2)=du1(3,jp,k)-du1(3,j,k) df3(3)=du1(3,j,k)-du1o(3,j,k) call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3) else jm=j-1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do df1(1)=du1(1,j,k)-du1(1,jm,k) df1(2)=du1(1,j,kp)-du1(1,j,km) df1(3)=du1(1,j,k)-du1o(1,j,k) df2(1)=du1(2,j,k)-du1(2,jm,k) df2(2)=du1(2,j,kp)-du1(2,j,km) df2(3)=du1(2,j,k)-du1o(2,j,k) df3(1)=du1(3,j,k)-du1(3,jm,k) df3(2)=du1(3,j,kp)-du1(3,j,km) df3(3)=du1(3,j,k)-du1o(3,j,k) call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3) end if c c derivatives of u c ux=du1(1,j,k) uy=du1(2,j,k) uz=du1(3,j,k) uxx=dgg1(1) uyy=dgg2(2) uzz=dgg3(3) uxy=(dgg1(2)+dgg2(1))/2.0d0 uxz=(dgg1(3)+dgg3(1))/2.0d0 uyz=(dgg2(3)+dgg3(2))/2.0d0 c c derivatives of S_I and Grad(S_I) c dfu=dffiu(j) dfuu=ddffiu(j) gx=ux*dfu gy=uy*dfu gz=uz*dfu gxx=dfuu*ux*ux+dfu*uxx gyy=dfuu*uy*uy+dfu*uyy gzz=dfuu*uz*uz+dfu*uzz gxy=dfuu*ux*uy+dfu*uxy gxz=dfuu*ux*uz+dfu*uxz gyz=dfuu*uy*uz+dfu*uyz c ggri(1,1,j,k)=gxx ggri(2,2,j,k)=gyy ggri(3,3,j,k)=gzz ggri(1,2,j,k)=gxy ggri(2,1,j,k)=gxy ggri(1,3,j,k)=gxz ggri(3,1,j,k)=gxz ggri(2,3,j,k)=gyz ggri(3,2,j,k)=gyz c c derivatives of |Grad(S_I)|^2 c dgrad2v(1,j,k)=2.0d0*(gx*gxx+gy*gxy+gz*gxz) dgrad2v(2,j,k)=2.0d0*(gx*gxy+gy*gyy+gz*gyz) dgrad2v(3,j,k)=2.0d0*(gx*gxz+gy*gyz+gz*gzz) end do end do c return end c c solution of the linear system of 3 eqs : dgg . dxv = dff c input vectors : dxv1, dxv2, dxv3, dff c output vector : dgg c c dff=(1,0,0) c subroutine solg0(dxv1,dxv2,dxv3,dgg) double precision denom,aa1,aa2,aa3 double precision dxv1(3),dxv2(3),dxv3(3),dgg(3) aa1=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) aa2=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) aa3=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) denom = dxv1(1)*aa1-dxv2(1)*aa2+dxv3(1)*aa3 dgg(1) = aa1/denom dgg(2) = -(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))/denom dgg(3) = (dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))/denom return end c c three rhs vectors df1, df2, df3 c subroutine solg3(dxv1,dxv2,dxv3,df1,df2,df3,dg1,dg2,dg3) double precision denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 double precision dxv1(3),dxv2(3),dxv3(3) double precision df1(3),df2(3),df3(3) double precision dg1(3),dg2(3),dg3(3) a11=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) a21=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) a31=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) a12=(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3)) a22=(dxv1(1)*dxv3(3)-dxv1(3)*dxv3(1)) a32=(dxv1(1)*dxv2(3)-dxv1(3)*dxv2(1)) a13=(dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2)) a23=(dxv1(1)*dxv3(2)-dxv1(2)*dxv3(1)) a33=(dxv1(1)*dxv2(2)-dxv1(2)*dxv2(1)) denom = dxv1(1)*a11-dxv2(1)*a21+dxv3(1)*a31 dg1(1) = (df1(1)*a11-df1(2)*a21+df1(3)*a31)/denom dg1(2) = -(df1(1)*a12-df1(2)*a22+df1(3)*a32)/denom dg1(3) = (df1(1)*a13-df1(2)*a23+df1(3)*a33)/denom dg2(1) = (df2(1)*a11-df2(2)*a21+df2(3)*a31)/denom dg2(2) = -(df2(1)*a12-df2(2)*a22+df2(3)*a32)/denom dg2(3) = (df2(1)*a13-df2(2)*a23+df2(3)*a33)/denom dg3(1) = (df3(1)*a11-df3(2)*a21+df3(3)*a31)/denom dg3(2) = -(df3(1)*a12-df3(2)*a22+df3(3)*a32)/denom dg3(3) = (df3(1)*a13-df3(2)*a23+df3(3)*a33)/denom return end c c Runge-Kutta integrator c subroutine rkint4 use beamdata, only : nrayr,nrayth,ywrk,ypwrk,grad2,dgrad2v, . gri,ggri implicit real*8 (a-h,o-z) parameter(ndim=6) dimension y(ndim),yy(ndim) dimension fk1(ndim),fk2(ndim),fk3(ndim),fk4(ndim) dimension dgr2(3),dgr(3),ddgr(3,3) c common/dsds/dst c common/gr/gr2 common/dgr/dgr2,dgr,ddgr common/igrad/igrad c h=dst hh=h*0.5d0 h6=h/6.0d0 c do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk do ieq=1,ndim y(ieq)=ywrk(ieq,j,k) fk1(ieq)=ypwrk(ieq,j,k) yy(ieq)=y(ieq)+fk1(ieq)*hh end do gr2=grad2(j,k) do iv=1,3 dgr2(iv)=dgrad2v(iv,j,k) dgr(iv)=gri(iv,j,k) do jv=1,3 ddgr(iv,jv)=ggri(iv,jv,j,k) end do end do call fwork(yy,fk2) c do ieq=1,ndim yy(ieq)=y(ieq)+fk2(ieq)*hh end do call fwork(yy,fk3) c do ieq=1,ndim yy(ieq)=y(ieq)+fk3(ieq)*h end do call fwork(yy,fk4) c do ieq=1,ndim ywrk(ieq,j,k)=y(ieq) . +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq)) end do end do end do c call updatepos c if(igrad.eq.1) call gradi c return end c c c subroutine gwork(j,k) use beamdata, only : ywrk,ypwrk,grad2,dgrad2v,gri,ggri implicit real*8 (a-h,o-z) parameter(ndim=6) dimension yy(ndim),yyp(ndim) dimension dgr2(3),dgr(3),ddgr(3,3) c common/igrad/igrad c common/gr/gr2 common/dgr/dgr2,dgr,ddgr c c begin --- update vector yy c do ieq=1,ndim yy(ieq)=ywrk(ieq,j,k) end do c if(igrad.eq.1) then gr2=grad2(j,k) do iv=1,3 dgr2(iv)=dgrad2v(iv,j,k) dgr(iv)=gri(iv,j,k) do jv=1,3 ddgr(iv,jv)=ggri(iv,jv,j,k) end do end do end if c call fwork(yy,yyp) c do ieq=1,ndim ypwrk(ieq,j,k)=yyp(ieq) end do c c end --- update vector yy c return end c c c subroutine fwork(y,dery) implicit real*8 (a-h,o-z) dimension y(6),dery(6) dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3) dimension derdxv(3),danpldxv(3),derdnv(3) dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3) c common/gr/gr2 common/dgr/dgr2,dgr,ddgr common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/mode/sox common/nplr/anpl,anpr common/bb/bv common/dbb/derbv common/igrad/igrad common/xgxg/xg common/ygyg/yg common/dxgyg/derxg,deryg common/vgv/vgm,derdnm common/dersdst/dersdst common/ierr/ierr common/anv/anv common/xv/xv common/idst/idst c xx=y(1) yy=y(2) zz=y(3) xv(1)=y(1) xv(2)=y(2) xv(3)=y(3) c anv(1) = y(4) anv(2) = y(5) anv(3) = y(6) c C rr=sqrt(xx**2+yy**2) C phi=acos(xx/rr) C if (yy.lt.0.0d0) then C phi=-phi C end if c call plas_deriv(xx,yy,zz) c an2=anv(1)*anv(1)+anv(2)*anv(2)+anv(3)*anv(3) anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3) c if(abs(anpl).gt.0.99d0) then if(abs(anpl).le.1.05d0) then ierr=97 else ierr=98 end if end if c anpl2=anpl*anpl dnl=1.0d0-anpl2 anpr2=an2-anpl2 anpr=0.0d0 if(anpr2.gt.0.0d0) anpr=sqrt(anpr2) yg2=yg**2 c an2s=1.0d0 dan2sdxg=0.0d0 dan2sdyg=0.0d0 dan2sdnpl=0.0d0 del=0.0d0 fdia=0.0d0 dfdiadnpl=0.0d0 dfdiadxg=0.0d0 dfdiadyg=0.0d0 c duh=1.0d0-xg-yg2 if(xg.gt.0.0d0) then del=sqrt(dnl*dnl+4.0d0*anpl*anpl*(1.0d0-xg)/yg2) an2s=1.0d0-xg -xg*yg2*(1.0d0+anpl2+sox*del)/duh/2.0d0 c dan2sdxg=-yg2*(1.0d0-yg2)*(1.0d0+anpl2+sox*del)/duh**2/2.0d0 . +sox*xg*anpl2/(del*duh)-1.0d0 dan2sdyg=-xg*yg*(1.0d0-xg)*(1.0d0+anpl2+sox*del)/duh**2 . +2.0d0*sox*xg*(1.0d0-xg)*anpl2/(yg*del*duh) dan2sdnpl=-xg*yg2*anpl/duh . -sox*xg*anpl*(2.0d0*(1.0d0-xg)-yg2*dnl)/(del*duh) c if(igrad.gt.0) then ddelnpl2=2.0d0*(2.0d0*(1.0d0-xg)*(1.0d0+3.0d0*anpl2**2) . -yg2*dnl**3)/yg2/del**3 fdia=-xg*yg2*(1.0d0+sox*ddelnpl2/2.0d0)/duh derdel=2.0d0*(1.0d0-xg)*anpl2*(1.0d0+3.0d0*anpl2**2) . -dnl**2*(1.0d0+3.0d0*anpl2)*yg2 derdel=4.0d0*derdel/(yg**5*del**5) ddelnpl2y=2.0d0*(1.0d0-xg)*derdel ddelnpl2x=yg*derdel dfdiadnpl=24.0d0*sox*xg*(1.0d0-xg)*anpl*(1.0d0-anpl**4) . /(yg2*del**5) dfdiadxg=-yg2*(1.0d0-yg2)/duh**2-sox*yg2*((1.0d0-yg2)*ddelnpl2 . +xg*duh*ddelnpl2x)/(2.0d0*duh**2) dfdiadyg=-2.0d0*yg*xg*(1.0d0-xg)/duh**2 . -sox*xg*yg*(2.0d0*(1.0d0-xg)*ddelnpl2 . +yg*duh*ddelnpl2y)/(2.0d0*duh**2) c end if end if c bdotgr=0.0d0 do iv=1,3 bdotgr=bdotgr+bv(iv)*dgr(iv) dbgr(iv)=0.0d0 do jv=1,3 dbgr(iv)=dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv) end do end do c derdnm=0.0d0 c do iv=1,3 danpldxv(iv)=anv(1)*derbv(1,iv)+anv(2)*derbv(2,iv) . +anv(3)*derbv(3,iv) derdxv(iv)=-(derxg(iv)*dan2sdxg . +deryg(iv)*dan2sdyg+danpldxv(iv)*dan2sdnpl) derdxv(iv)=derdxv(iv)-igrad*dgr2(iv) c derdxv(iv)=derdxv(iv)+fdia*bdotgr*dbgr(iv)+0.5d0*bdotgr**2* . (derxg(iv)*dfdiadxg+deryg(iv)*dfdiadyg+danpldxv(iv)*dfdiadnpl) c derdnv(iv)=2.0d0*anv(iv)-bv(iv)*dan2sdnpl . +0.5d0*bdotgr**2*bv(iv)*dfdiadnpl c derdnm=derdnm+derdnv(iv)**2 end do c derdnm=sqrt(derdnm) c derdom=-2.0d0*an2+2.0d0*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl . +2.0d0*igrad*gr2 . -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0d0 . +anpl*dfdiadnpl/2.0d0) c if (idst.eq.0) then c integration variable: s denom=-derdnm else if (idst.eq.1) then c integration variable: c*t denom=derdom else c integration variable: Sr denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3)) end if c c coefficient for integration in s c ds/dst, where st is the integration variable dersdst=-derdnm/denom c dery(1) = -derdnv(1)/denom dery(2) = -derdnv(2)/denom dery(3) = -derdnv(3)/denom dery(4) = derdxv(1)/denom dery(5) = derdxv(2)/denom dery(6) = derdxv(3)/denom c c vgv : ~ group velocity c vgm=0 do iv=1,3 vgv(iv)=-derdnv(iv)/derdom vgm=vgm+vgv(iv)**2 end do vgm=sqrt(vgm) c c dd : dispersion relation (real part) c ddi : dispersion relation (imaginary part) c dd=an2-an2s-igrad*(gr2-0.5d0*bdotgr**2*fdia) ddi=derdnv(1)*dgr(1)+derdnv(2)*dgr(2)+derdnv(3)*dgr(3) c return end c c c subroutine plas_deriv(xx,yy,zz) use const_and_precisions, only : pi implicit real*8 (a-h,o-z) dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3) dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3) c common/parbres/bres common/parpl/brr,bphi,bzz,ajphi common/btot/btot common/bb/bv common/dbb/derbv common/xgxg/xg common/dxgdps/dxgdpsi common/ygyg/yg common/dxgyg/derxg,deryg common/iieq/iequil common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/psinv/psinv common/sgnib/sgnbphi,sgniphi c xg=0.0d0 yg=9.9d1 c do iv=1,3 derxg(iv)=0.0d0 deryg(iv)=0.0d0 bv(iv)=0.0d0 dbtot(iv)=0.0d0 do jv=1,3 dbv(iv,jv)=0.0d0 derbv(iv,jv)=0.0d0 dbvcdc(iv,jv)=0.0d0 dbvcdc(iv,jv)=0.0d0 dbvdc(iv,jv)=0.0d0 end do end do c if(iequil.eq.0) return c c cylindrical coordinates c rr2=xx**2+yy**2 rr=sqrt(rr2) csphi=xx/rr snphi=yy/rr c bv(1)=-snphi*sgnbphi bv(2)=csphi*sgnbphi c c convert from cm to meters c zzm=1.0d-2*zz rrm=1.0d-2*rr c if(iequil.eq.1) then call equian(rrm,zzm) yg=fpolv/(rrm*bres) end if c if(iequil.eq.2) then call equinum(rrm,zzm) call sub_xg_derxg yg=fpolv/(rrm*bres) bphi=fpolv/rrm btot=abs(bphi) if(psinv.lt.0.0d0) return end if c c B = f(psi)/R e_phi+ grad psi x e_phi/R c c bvc(i) = B_i in cylindrical coordinates c bv(i) = B_i in cartesian coordinates c bphi=fpolv/rrm brr=-dpsidz/rrm bzz= dpsidr/rrm c dfpolv=ffpv/fpolv c bvc(1)=brr bvc(2)=bphi bvc(3)=bzz c bv(1)=bvc(1)*csphi-bvc(2)*snphi bv(2)=bvc(1)*snphi+bvc(2)*csphi bv(3)=bvc(3) c b2tot=bv(1)**2+bv(2)**2+bv(3)**2 btot=sqrt(b2tot) c c dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) c dbvcdc(1,1)=-ddpsidrz/rrm-brr/rrm dbvcdc(1,3)=-ddpsidzz/rrm dbvcdc(2,1)= dfpolv*dpsidr/rrm-bphi/rrm dbvcdc(2,3)= dfpolv*dpsidz/rrm dbvcdc(3,1)= ddpsidrr/rrm-bzz/rrm dbvcdc(3,3)= ddpsidrz/rrm c c dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) c dbvdc(1,1) = dbvcdc(1,1)*csphi-dbvcdc(2,1)*snphi dbvdc(1,2) = -bv(2) dbvdc(1,3) = dbvcdc(1,3)*csphi-dbvcdc(2,3)*snphi dbvdc(2,1) = dbvcdc(1,1)*snphi+dbvcdc(2,1)*csphi dbvdc(2,2) = bv(1) dbvdc(2,3) = dbvcdc(1,3)*snphi+dbvcdc(2,3)*csphi dbvdc(3,1) = dbvcdc(3,1) dbvdc(3,2) = dbvcdc(3,2) dbvdc(3,3) = dbvcdc(3,3) c drrdx=csphi drrdy=snphi dphidx=-snphi/rrm dphidy=csphi/rrm c c dbv(iv,jv) = d Bcart(iv) / dxvcart(jv) c dbv(1,1)=drrdx*dbvdc(1,1)+dphidx*dbvdc(1,2) dbv(1,2)=drrdy*dbvdc(1,1)+dphidy*dbvdc(1,2) dbv(1,3)=dbvdc(1,3) dbv(2,1)=drrdx*dbvdc(2,1)+dphidx*dbvdc(2,2) dbv(2,2)=drrdy*dbvdc(2,1)+dphidy*dbvdc(2,2) dbv(2,3)=dbvdc(2,3) dbv(3,1)=drrdx*dbvdc(3,1)+dphidx*dbvdc(3,2) dbv(3,2)=drrdy*dbvdc(3,1)+dphidy*dbvdc(3,2) dbv(3,3)=dbvdc(3,3) c dbtot(1)=(bv(1)*dbv(1,1)+bv(2)*dbv(2,1)+bv(3)*dbv(3,1))/btot dbtot(2)=(bv(1)*dbv(1,2)+bv(2)*dbv(2,2)+bv(3)*dbv(3,2))/btot dbtot(3)=(bv(1)*dbv(1,3)+bv(2)*dbv(2,3)+bv(3)*dbv(3,3))/btot c yg=btot/Bres c c convert spatial derivatives from dummy/m -> dummy/cm c to be used in fwork c c bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z c do iv=1,3 deryg(iv)=1.0d-2*dbtot(iv)/Bres bv(iv)=bv(iv)/btot do jv=1,3 derbv(iv,jv)=1.0d-2*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot end do end do c derxg(1)=1.0d-2*drrdx*dpsidr*dxgdpsi derxg(2)=1.0d-2*drrdy*dpsidr*dxgdpsi derxg(3)=1.0d-2*dpsidz*dxgdpsi c c current density computation in Ampere/m^2 c ccj==1/mu_0 c ccj=1.0d+7/(4.0d0*pi) ajphi=ccj*(dbvcdc(1,3)-dbvcdc(3,1)) c ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3)) c ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2)) c return end c c c subroutine equian(rrm,zzm) implicit real*8 (a-h,o-z) c common/parqq/q0,qa,alq common/parban/b0,rr0m,zr0m,rpam common/psinv/psinv common/pareq1/psia common/densbnd/psdbnd common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn common/dddens/dens,ddens common/sgnib/sgnbphi,sgniphi common/bmxmn/bmxi,bmni common/fc/fci common/iipr/iprof c if(iprof.eq.0) psdbnd=1.0d0 c dpsidrp=0.0d0 d2psidrp=0.0d0 c c simple model for equilibrium c rpm=sqrt((rrm-rr0m)**2+(zzm-zr0m)**2) rn=rpm/rpam c snt=0.0d0 cst=1.0d0 if (rpm.gt.0.0d0) then snt=(zzm-zr0m)/rpm cst=(rrm-rr0m)/rpm end if c c valore fittizio di psinv e di psia: psinv=rn**2 psia=sgniphi c sgn=sgniphi*sgnbphi if(rn.le.1.0d0) then qq=q0+(qa-q0)*rn**alq dpsidrp=B0*rpam*rn/qq*sgn dqq=alq*(qa-q0)*rn**(alq-1.0d0) d2psidrp=B0*(1.0d0-rn*dqq/qq)/qq*sgn else dpsidrp=B0*rpam/qa/rn*sgn d2psidrp=-B0/qa/rn**2 *sgn end if c fpolv=sgnbphi*b0*rr0m dfpolv=0.0d0 ffpv=sgniphi*fpolv*dfpolv c dpsidr=dpsidrp*cst dpsidz=dpsidrp*snt ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm) ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp c dens=0.0d0 xg=0.0d0 dxgdpsi=0.0d0 if(psinv.lt.psdbnd) then call density(psinv) end if xg=xgcn*dens ddensdrp=2.0d0*rn*ddens if(dpsidrp.ne.0.0d0) dxgdpsi=xgcn*ddensdrp/(dpsidrp*rpam) c bmax=sqrt(fpolv**2+dpsidrp**2)/(rr0m-rpam*rn) bmin=sqrt(fpolv**2+dpsidrp**2)/(rr0m+rpam*rn) c return end c c c subroutine equinum(rpsim,zpsim) implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) parameter(lwrk=8,liwrk=2) dimension ccspl(nnw*nnh) dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) dimension rrs(1),zzs(1),ffspl(1) parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) parameter(lw11=nnw*3+nnh*3+nnw*nnh) parameter(nrs=1,nzs=1) dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) dimension tfp(nrest),cfp(nrest),wrkfd(nrest) c common/eqnn/nr,nz,nrho,npp,nintp common/psinv/psinv common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/pareq1/psia common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/pareq1t/psiant,psinop c common/coffeqt/tr,tz common/coffeqtp/tfp common/coffeq/ccspl common/coffeqd/cc01,cc10,cc20,cc02,cc11 common/coffeqn/nsrt,nszt,nsft common/cofffp/cfp common/fpas/fpolas c psinv=-1.0d0 c dpsidr=0.0d0 dpsidz=0.0d0 ddpsidrr=0.0d0 ddpsidzz=0.0d0 ddpsidrz=0.0d0 c c here lengths are measured in meters c fpolv=fpolas ffpv=0.0d0 c if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return c rrs(1)=rpsim zzs(1)=zpsim nsr=nsrt nsz=nszt call dierckx_fpbisp(tr,nsr,tz,nsz,ccspl,3,3, . rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2)) psinv=(ffspl(1)-psinop)/psiant if(psinv.lt.0.0d0) . write(*,'(a,3e12.4)') ' psin < 0 , R, z ',psinv,rpsim,zpsim c nur=1 nuz=0 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10, . kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc10(iwr),cc10(iwz),iwrk(1),iwrk(2)) dpsidr= ffspl(1)*psia c nur=0 nuz=1 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01, . kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc01(iwr),cc01(iwz),iwrk(1),iwrk(2)) dpsidz= ffspl(1)*psia c nur=2 nuz=0 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20, . kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc20(iwr),cc20(iwz),iwrk(1),iwrk(2)) ddpsidrr= ffspl(1)*psia c nur=0 nuz=2 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02, . kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc02(iwr),cc02(iwz),iwrk(1),iwrk(2)) ddpsidzz= ffspl(1)*psia c nur=1 nuz=1 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11, . kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc11(iwr),cc11(iwz),iwrk(1),iwrk(2)) ddpsidrz= ffspl(1)*psia c if(psinv.le.1.0d0.and.psinv.gt.0.0d0) then rrs(1)=psinv call dierckx_splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier) fpolv=ffspl(1) nu=1 call dierckx_splder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier) dfpolv=ffspl(1) ffpv=fpolv*dfpolv/psia end if c return end c c c subroutine bfield(rpsim,zpsim,bphi,brr,bzz) implicit real*8 (a-h,o-z) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv call equinum(rpsim,zpsim) bphi=fpolv/rpsim brr=-dpsidz/rpsim bzz= dpsidr/rpsim return end c subroutine tor_curr_psi(h,ajphi) implicit real*8 (a-h,o-z) common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz call psi_raxis(h,r1,r2) call tor_curr(r2,zmaxis,ajphi) return end c subroutine tor_curr(rpsim,zpsim,ajphi) use const_and_precisions, only : pi, mu0=>mu0_ implicit real*8 (a-h,o-z) parameter(ccj=1.0d0/mu0) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv call equinum(rpsim,zpsim) bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim ajphi=ccj*(dbvcdc13-dbvcdc31) return end c subroutine psi_raxis(h,r1,r2) implicit real*8 (a-h,o-z) parameter(mest=4,kspl=3) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) dimension cc(nnw*nnh),tr(nrest),tz(nzest) dimension czc(nrest),zeroc(mest) c common/pareq1t/psiant,psinop common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/coffeqn/nsr,nsz,nsft common/coffeq/cc common/coffeqt/tr,tz c iopt=1 zc=zmaxis call dierckx_sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc, . ier) if(ier.gt.0) write(*,*) ' dierckx_sprofil =',ier val=h*psiant+psinop call dierckx_sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) r2=zeroc(2) return end c c c subroutine sub_xg_derxg implicit real*8 (a-h,o-z) common/psinv/psinv common/pareq1/psia common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn common/dddens/dens,ddenspsin xg=0.0d0 dxgdpsi=0.0d0 c if(psinv.le.psdbnd.and.psinv.ge.0) then call density(psinv) xg=xgcn*dens dxgdpsi=xgcn*ddenspsin/psia c end if return end c c c subroutine density(arg) implicit real*8 (a-h,o-z) parameter(npmx=501,npest=npmx+4) dimension xxs(1),ffs(1) dimension tfn(npest),cfn(npest),wrkfd(npest) c common/densbnd/psdbnd common/denspp/psnpp,aad,bbd,ccd common/iipr/iprof common/pardens/dens0,aln1,aln2 common/dddens/dens,ddens common/coffdt/tfn common/coffdnst/nsfd common/cofffn/cfn c c computation of density [10^19 m^-3] and derivative wrt psi c dens=0.0d0 ddens=0.0d0 if(arg.gt.psdbnd.or.arg.lt.0.0d0) return c if(iprof.eq.0) then if(arg.gt.1.0d0) return profd=(1.0d0-arg**aln1)**aln2 dens=dens0*profd dprofd=-aln1*aln2*arg**(aln1-1.0d0) . *(1.0d0-arg**aln1)**(aln2-1.0d0) ddens=dens0*dprofd else if(arg.le.psdbnd.and.arg.gt.psnpp) then c c cubic interpolation for 1 < psi < psdbnd c nn=3 nn1=nn+1 nn2=nn+2 dpsib=arg-psdbnd dens=aad*dpsib**nn+bbd*dpsib**nn1+ccd*dpsib**nn2 ddens=nn*aad*dpsib**(nn-1)+nn1*bbd*dpsib**nn . +nn2*ccd*dpsib**nn1 else xxs(1)=arg ier=0 call dierckx_splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier) dens=ffs(1) nu=1 ier=0 call dierckx_splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier) ddens=ffs(1) if(ier.gt.0) write(*,*) ier if(abs(dens).lt.1.0d-10) dens=0.0d0 end if if(dens.lt.0.0d0) write(*,*) ' DENSITY NEGATIVE',dens c end if return end c c c function temperature(arg) implicit real*8 (a-h,o-z) parameter(npmx=501) dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4) c common/parqte/te0,dte0,alt1,alt2 common/iipr/iprof common/crad/psrad,derad,terad,zfc common/coffte/ct common/eqnn/nr,nz,nrho,npp,nintp c temperature=0.0d0 if(arg.ge.1.0d0.or.arg.lt.0.0d0) return if(iprof.eq.0) then proft=(1.0d0-arg**alt1)**alt2 temperature=(te0-dte0)*proft+dte0 else call vlocate(psrad,npp,arg,k) k=max(1,min(k,npp-1)) dps=arg-psrad(k) temperature=fspli(ct,npmx,k,dps) endif return end c c c function fzeff(arg) implicit real*8 (a-h,o-z) parameter(npmx=501) dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),cz(npmx,4) c common/iipr/iprof common/crad/psrad,derad,terad,zfc common/coffz/cz common/eqnn/nr,nz,nrho,npp,nintp common/zz/Zeff c fzeff=1 ps=arg if(ps.gt.1.0d0.and.ps.lt.0.0d0) return if(iprof.eq.0) then fzeff=zeff else call vlocate(psrad,npp,ps,k) k=max(1,min(k,npp-1)) dps=ps-psrad(k) fzeff=fspli(cz,npmx,k,dps) endif return end c c beam tracing initial conditions igrad=1 c subroutine ic_gb use const_and_precisions, only : pi use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v, . gri,ggri implicit real*8 (a-h,o-z) parameter(zero=0.0d0,izero=0,one=1.0d0) parameter(cvdr=pi/180.0d0) dimension ytmp(6),yptmp(6) complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy complex*16 catand external catand c parameter(ui=(0.0d0,1.0d0)) c common/rwmax/rwmax common/parwv/ak0,akinv,fhz common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/psinv/psinv common/dddens/dens,ddens common/tete/tekev common/nplr/anpl,anpr common/parpl/brr,bphi,bzz,ajphi c ui=(0.0d0,1.0d0) csth=anz0c snth=sqrt(1.0d0-csth**2) csps=1.0d0 snps=0.0d0 if(snth.gt.0.0d0) then csps=any0c/snth snps=anx0c/snth end if c phiwrad=phiw*cvdr phirrad=phir*cvdr csphiw=cos(phiwrad) snphiw=sin(phiwrad) c csphir=cos(phirrad) c snphir=sin(phirrad) c wwcsi=2.0d0*akinv/wcsi**2 wweta=2.0d0*akinv/weta**2 c if(phir.ne.phiw) then sk=(rcicsi+rcieta) sw=(wwcsi+wweta) dk=(rcicsi-rcieta) dw=(wwcsi-wweta) ts=-(dk*sin(2.0d0*phirrad)-ui*dw*sin(2.0d0*phiwrad)) tc=(dk*cos(2.0d0*phirrad)-ui*dw*cos(2.0d0*phiwrad)) phic=0.5d0*catand(ts/tc) ddd=dk*cos(2*(phirrad+phic))-ui*dw*cos(2*(phiwrad+phic)) sss=sk-ui*sw qi1=0.5d0*(sss+ddd) qi2=0.5d0*(sss-ddd) rci1=dble(qi1) ww1=-dimag(qi1) rci2=dble(qi2) ww2=-dimag(qi2) else rci1=rcicsi rci2=rcieta ww1=wwcsi ww2=wweta phic=-phiwrad qi1=rci1-ui*ww1 qi2=rci2-ui*ww2 end if c w01=sqrt(2.0d0*akinv/ww1) c z01=-rci1/(rci1**2+ww1**2) c w02=sqrt(2.0d0*akinv/ww2) c z02=-rci2/(rci2**2+ww2**2) qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2 qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2 qqxy=-(qi1-qi2)*sin(2.0d0*phic) wwxx=-dimag(qqxx) wwyy=-dimag(qqyy) wwxy=-dimag(qqxy)/2.0d0 rcixx=dble(qqxx) rciyy=dble(qqyy) rcixy=dble(qqxy)/2.0d0 dqi1=-qi1**2 dqi2=-qi2**2 d2qi1=2*qi1**3 d2qi2=2*qi2**3 dqqxx=dqi1*cos(phic)**2+dqi2*sin(phic)**2 dqqyy=dqi1*sin(phic)**2+dqi2*cos(phic)**2 dqqxy=-(dqi1-dqi2)*sin(2.0d0*phic) d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2 d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2 d2qqxy=-(d2qi1-d2qi2)*sin(2.0d0*phic) dwwxx=-dimag(dqqxx) dwwyy=-dimag(dqqyy) dwwxy=-dimag(dqqxy)/2.0d0 d2wwxx=-dimag(d2qqxx) d2wwyy=-dimag(d2qqyy) d2wwxy=-dimag(d2qqxy)/2.0d0 drcixx=dble(dqqxx) drciyy=dble(dqqyy) drcixy=dble(dqqxy)/2.0d0 dr=1.0d0 if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) da=2.0d0*pi/dble(nrayth) c ddfu=2.0d0*dr**2*akinv do j=1,nrayr u=dble(j-1) c ffi=u**2*ddfu/2.0d0 dffiu(j)=u*ddfu ddffiu(j)=ddfu do k=1,nrayth alfak=(k-1)*da dcsiw=dr*cos(alfak)*wcsi detaw=dr*sin(alfak)*weta dx0t=dcsiw*csphiw-detaw*snphiw dy0t=dcsiw*snphiw+detaw*csphiw x0t=u*dx0t y0t=u*dy0t z0t=-0.5d0*(rcixx*x0t**2+rciyy*y0t**2+2*rcixy*x0t*y0t) c c csiw=u*dcsiw c etaw=u*detaw c csir=x0t*csphir+y0t*snphir c etar=-x0t*snphir+y0t*csphir c dx0= x0t*csps+snps*(y0t*csth+z0t*snth) dy0=-x0t*snps+csps*(y0t*csth+z0t*snth) dz0= z0t*csth-y0t*snth x0=x00+dx0 y0=y00+dy0 z0=z00+dz0 gxt=x0t*wwxx+y0t*wwxy gyt=x0t*wwxy+y0t*wwyy gzt=0.5d0*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy gr2=gxt*gxt+gyt*gyt+gzt*gzt gxxt=wwxx gyyt=wwyy gzzt=0.5d0*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy gxyt=wwxy gxzt=x0t*dwwxx+y0t*dwwxy gyzt=x0t*dwwxy+y0t*dwwyy dgr2xt=2.0d0*(gxt*gxxt+gyt*gxyt+gzt*gxzt) dgr2yt=2.0d0*(gxt*gxyt+gyt*gyyt+gzt*gyzt) dgr2zt=2.0d0*(gxt*gxzt+gyt*gyzt+gzt*gzzt) dgr2x= dgr2xt*csps+snps*(dgr2yt*csth+dgr2zt*snth) dgr2y=-dgr2xt*snps+csps*(dgr2yt*csth+dgr2zt*snth) dgr2z= dgr2zt*csth-dgr2yt*snth gri(1,j,k)=gxt*csps+snps*(gyt*csth+gzt*snth) gri(2,j,k)=-gxt*snps+csps*(gyt*csth+gzt*snth) gri(3,j,k)=gzt*csth-gyt*snth ggri(1,1,j,k)=gxxt*csps**2+ . snps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)+ . 2.0d0*snps*csps*(gxyt*csth+gxzt*snth) ggri(2,2,j,k)=gxxt*snps**2+ . csps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)- . 2.0d0*snps*csps*(gxyt*csth+gxzt*snth) ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2-2.0d0*csth*snth*gyzt ggri(1,2,j,k)=csps*snps*(-gxxt+csth**2*gyyt+snth**2*gzzt . +2.0d0*csth*snth*gyzt) . +(csps**2-snps**2)*(snth*gxzt+csth*gxyt) ggri(1,3,j,k)=csth*snth*snps*(gzzt-gyyt) . +(csth**2-snth**2)*snps*gyzt+csps*(csth*gxzt-snth*gxyt) ggri(2,3,j,k)=csth*snth*csps*(gzzt-gyyt) . +(csth**2-snth**2)*csps*gyzt+snps*(snth*gxyt-csth*gxzt) ggri(2,1,j,k)=ggri(1,2,j,k) ggri(3,1,j,k)=ggri(1,3,j,k) ggri(3,2,j,k)=ggri(2,3,j,k) c du1tx=(dx0t*wwxx+dy0t*wwxy)/ddfu du1ty=(dx0t*wwxy+dy0t*wwyy)/ddfu du1tz=0.5d0*u*(dx0t**2*dwwxx+dy0t**2*dwwyy+ . 2.0d0*dx0t*dy0t*dwwxy)/ddfu c pppx=x0t*rcixx+y0t*rcixy pppy=x0t*rcixy+y0t*rciyy denpp=pppx*gxt+pppy*gyt if (denpp.ne.0.0d0) then ppx=-pppx*gzt/denpp ppy=-pppy*gzt/denpp else ppx=0.0d0 ppy=0.0d0 end if c anzt=sqrt((1.0d0+gr2)/(1.0d0+ppx**2+ppy**2)) anxt=ppx*anzt anyt=ppy*anzt c anx= anxt*csps+snps*(anyt*csth+anzt*snth) any=-anxt*snps+csps*(anyt*csth+anzt*snth) anz= anzt*csth-anyt*snth c an20=1.0d0+gr2 an0=sqrt(an20) anx0=anx any0=any anz0=anz c xc0(1,j,k)=x0 xc0(2,j,k)=y0 xc0(3,j,k)=z0 c ywrk0(1,j,k)=x0 ywrk0(2,j,k)=y0 ywrk0(3,j,k)=z0 ywrk0(4,j,k)=anx0 ywrk0(5,j,k)=any0 ywrk0(6,j,k)=anz0 c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 ypwrk0(4,j,k) = dgr2x/an0/2.0d0 ypwrk0(5,j,k) = dgr2y/an0/2.0d0 ypwrk0(6,j,k) = dgr2z/an0/2.0d0 c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) grad2(j,k)=gr2 dgrad2v(1,j,k)=dgr2x dgrad2v(2,j,k)=dgr2y dgrad2v(3,j,k)=dgr2z c du10(1,j,k)= du1tx*csps+snps*(du1ty*csth+du1tz*snth) du10(2,j,k)=-du1tx*snps+csps*(du1ty*csth+du1tz*snth) du10(3,j,k)= du1tz*csth-du1ty*snth c dd=anx0**2+any0**2+anz0**2-an20 vgradi=anxt*gxt+anyt*gyt+anzt*gzt ddi=2.0d0*vgradi c r0=sqrt(x0**2+y0**2) x0m=x0/1.0d2 y0m=y0/1.0d2 r0m=r0/1.0d2 z0m=z0/1.0d2 if(j.eq.nrayr) then write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, . psinv,one,dens,tekev,brr,bphi,bzz, . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one ddr110=dd end if if(j.eq.nrayr.and.k.eq.1) then write(617,99) zero,ddr110,dd,ddi end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=608 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end c c ray tracing initial conditions igrad=0 c subroutine ic_rt use const_and_precisions, only : pi use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v, . gri,ggri implicit real*8 (a-h,o-z) parameter(zero=0.0d0,izero=0,one=1.0d0) parameter(cvdr=pi/180.0d0) dimension ytmp(6),yptmp(6) c common/rwmax/rwmax common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/psinv/psinv common/dddens/dens,ddens common/tete/tekev common/nplr/anpl,anpr common/parpl/brr,bphi,bzz,ajphi c csth=anz0c snth=sqrt(1.0d0-csth**2) csps=1.0d0 snps=0.0d0 if(snth.gt.0.0d0) then csps=any0c/snth snps=anx0c/snth end if c phiwrad=phiw*cvdr csphiw=cos(phiwrad) snphiw=sin(phiwrad) c dr=1.0d0 if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) da=2.0d0*pi/dble(nrayth) z0t=0.0d0 c do j=1,nrayr u=dble(j-1) dffiu(j)=0.0d0 ddffiu(j)=0.0d0 do k=1,nrayth alfak=(k-1)*da dcsiw=dr*cos(alfak)*wcsi detaw=dr*sin(alfak)*weta dx0t=dcsiw*csphiw-detaw*snphiw dy0t=dcsiw*snphiw+detaw*csphiw x0t=u*dx0t y0t=u*dy0t c c csiw=u*dcsiw c etaw=u*detaw c csir=csiw c etar=etaw c dx0= x0t*csps+snps*(y0t*csth+z0t*snth) dy0=-x0t*snps+csps*(y0t*csth+z0t*snth) dz0= z0t*csth-y0t*snth c x0=x00+dx0 y0=y00+dy0 z0=z00+dz0 c ppcsi=u*dr*cos(alfak)*rcicsi ppeta=u*dr*sin(alfak)*rcieta c anzt=1.0d0/sqrt(1.0d0+ppcsi**2+ppeta**2) ancsi=ppcsi*anzt aneta=ppeta*anzt c anxt=ancsi*csphiw-aneta*snphiw anyt=ancsi*snphiw+aneta*csphiw c anx= anxt*csps+snps*(anyt*csth+anzt*snth) any=-anxt*snps+csps*(anyt*csth+anzt*snth) anz= anzt*csth-anyt*snth c an20=1.0d0 an0=sqrt(an20) anx0=anx any0=any anz0=anz c xc0(1,j,k)=x0 xc0(2,j,k)=y0 xc0(3,j,k)=z0 c ywrk0(1,j,k)=x0 ywrk0(2,j,k)=y0 ywrk0(3,j,k)=z0 ywrk0(4,j,k)=anx0 ywrk0(5,j,k)=any0 ywrk0(6,j,k)=anz0 c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) do iv=1,3 gri(iv,j,k)=0.0d0 dgrad2v(iv,j,k)=0.0d0 du10(iv,j,k)=0.0d0 do jv=1,3 ggri(iv,jv,j,k)=0.0d0 end do end do grad2(j,k)=0.0d0 c dd=anx0**2+any0**2+anz0**2-an20 vgradi=0.0d0 ddi=2.0d0*vgradi c r0=sqrt(x0**2+y0**2) x0m=x0/1.0d2 y0m=y0/1.0d2 r0m=r0/1.0d2 z0m=z0/1.0d2 if(j.eq.nrayr) then write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(617,99) zero,zero,zero,zero write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, . psinv,one,dens,tekev,brr,bphi,bzz, . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=608 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end subroutine ic_rt2 use const_and_precisions, only : pi use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,grad2,dgrad2v, . gri,ggri,yyrfl implicit real*8 (a-h,o-z) parameter(zero=0.0d0,izero=0,one=1.0d0) parameter(cvdr=pi/180.0d0) dimension ytmp(6),yptmp(6) c common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/psinv/psinv common/dddens/dens,ddens common/tete/tekev common/nplr/anpl,anpr common/parpl/brr,bphi,bzz,ajphi do j=1,nrayr do k=1,nrayth x0=yyrfl(j,k,1) y0=yyrfl(j,k,2) z0=yyrfl(j,k,3) anx0=yyrfl(j,k,4) any0=yyrfl(j,k,5) anz0=yyrfl(j,k,6) an20=anx0*anx0+any0*any0+anz0*anz0 an0=sqrt(an20) c xc0(1,j,k)=x0 xc0(2,j,k)=y0 xc0(3,j,k)=z0 c ywrk0(1,j,k)=x0 ywrk0(2,j,k)=y0 ywrk0(3,j,k)=z0 ywrk0(4,j,k)=anx0/an0 ywrk0(5,j,k)=any0/an0 ywrk0(6,j,k)=anz0/an0 c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) do iv=1,3 gri(iv,j,k)=0.0d0 dgrad2v(iv,j,k)=0.0d0 du10(iv,j,k)=0.0d0 do jv=1,3 ggri(iv,jv,j,k)=0.0d0 end do end do grad2(j,k)=0.0d0 c dd=anx0**2+any0**2+anz0**2-an20 vgradi=0.0d0 ddi=2.0d0*vgradi c r0=sqrt(x0**2+y0**2) x0m=x0/1.0d2 y0m=y0/1.0d2 r0m=r0/1.0d2 z0m=z0/1.0d2 if(j.eq.nrayr) then write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(617,99) zero,zero,zero,zero write(604,99) strfl11,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, . psinv,one,dens,tekev,brr,bphi,bzz, . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=608 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end c c ray power weigth coefficient q(j) c subroutine pweigth use beamdata, only : nrayr,nrayth,q implicit real*8(a-h,o-z) common/rwmax/rwmax c dr=1.0d0 if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) r1=0.0d0 summ=0.0d0 q(1)=1.0d0 if(nrayr.gt.1) then do j=1,nrayr r2=(dble(j)-0.5d0)*dr q(j)=(exp(-2.0d0*r1**2)-exp(-2.0d0*r2**2)) r1=r2 summ=summ+q(j) end do c q(1)=q(1)/summ sm=q(1) j=1 k=1 do j=2,nrayr q(j)=q(j)/nrayth/summ do k=1,nrayth sm=sm+q(j) end do end do end if c return end c c c subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi, . bmxi,bmni,fci,intp) implicit real*8 (a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4) dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4) dimension carea(nnintp,4),cfc(nnintp,4) c common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc common/pstab/rpstab common/eqnn/nr,nz,nrho,npp,nintp c ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) c dps=rpsi-rpstab(ip) c areai=fspli(carea,nnintp,ip,dps) voli=fspli(cvol,nnintp,ip,dps) dervoli=fsplid(cvol,nnintp,ip,dps) rrii=fspli(crri,nnintp,ip,dps) c if(intp.eq.0) return c rbavi=fspli(crbav,nnintp,ip,dps) bmxi=fspli(cbmx,nnintp,ip,dps) bmni=fspli(cbmn,nnintp,ip,dps) fci=fspli(cfc,nnintp,ip,dps) c return end c c c subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) implicit real*8 (a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp) dimension cratja(nnintp,4),cratjb(nnintp,4),cratjpl(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,nrho,npp,nintp common/cratj/cratja,cratjb,cratjpl ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) dps=rpsi-rpstab(ip) ratjai=fspli(cratja,nnintp,ip,dps) ratjbi=fspli(cratjb,nnintp,ip,dps) ratjpli=fspli(cratjpl,nnintp,ip,dps) return end c c c subroutine pabs_curr(i,j,k) use const_and_precisions, only : pi use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,didst, . ccci,q,tau1v implicit real*8 (a-h,o-z) c common/p0/p0mw c common/dsds/dst common/dersdst/dersdst common/iieq/iequil c common/parban/b0,rr0m,zr0m,rpam common/absor/alpha,effjcd,akim,tau0 c common/psinv/psinv common/sgnib/sgnbphi,sgniphi common/bmxmn/bmxi,bmni common/fc/fci c dvvl=1.0d0 rbavi=1.0d0 rrii=rr0m rhop0=sqrt(psjki(j,k,i-1)) if(iequil.eq.2) then intp=1 psinv=psjki(j,k,i) rhop=sqrt(psinv) call valpsispl(rhop,voli,dervoli,area,rrii, . rbavi,bmxi,bmni,fci,intp) dvvl=dervoli*abs(rhop-rhop0) else dvvl=4.0d0*pi*rr0m*pi*abs(rhop*rhop-rhop0*rhop0)*rpam**2 rbavi=b0/bmni fci=1.0d0-1.469d0*sqrt(rpam/rr0m) end if if(dvvl.le.0.0d0) dvvl=1.0d0 c adnm=2.0d0*pi*rrii c c calcolo della corrente cohen: currtot [MA] c calcolo della densita' corrente cohen: currj [MA/m^2] c calcolo della densita' potenza: pdjki [MW/m^3] c calcolo della efficienza : effjcdav [A m/W ] c tau0=tauv(j,k,i-1) c call ecrh_cd c alphav(j,k,i)=alpha dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0 tauv(j,k,i)=tauv(j,k,i-1)+dtau dpdst=p0mw*q(j)*exp(-tauv(j,k,i)-tau1v(j,k))* . alphav(j,k,i)*dersdst pdjki(j,k,i)=dpdst*dst/dvvl ppabs(j,k,i)=p0mw*q(j)*exp(-tau1v(j,k))* . (1.0d0-exp(-tauv(j,k,i))) effjcdav=rbavi*effjcd currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i) didst(j,k,i)=sgnbphi*effjcdav*dpdst/adnm dcidst=(didst(j,k,i)+didst(j,k,i-1))/2.0d0 ccci(j,k,i)=ccci(j,k,i-1)+dcidst*dst return end c c c subroutine ecrh_cd implicit real*8 (a-h,o-z) real*8 mc2,mesi parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8) parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) parameter(vcsi=2.99792458d+8) parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3) c common/ithn/ithn common/nharm/nharm,nhf common/warm/iwarm,ilarm common/ieccd/ieccd c common/ygyg/yg common/nplr/anpl,anpr common/vgv/vgm,derdnm common/parwv/ak0,akinv,fhz c common/nprw/anprre,anprim c common/absor/alpha,effjcd,akim,tau c common/psinv/psinv common/tete/tekev common/amut/amu common/zz/Zeff c c absorption computation c alpha=0.0d0 akim=0.0d0 effjcd=0.0d0 c tekev=temperature(psinv) amu=mc2/tekev zeff=fzeff(psinv) c if(tekev.le.0.0d0.or.tau.gt.taucr) return c dnl=1.0d0-anpl*anpl if(iwarm.eq.1) then ygc=1.0d0-anpl**2/2.0d0 else ygc=sqrt(1.0d0-anpl**2) end if c nharm=int(ygc/yg-eps)+1 c if(nharm.eq.0) return c nhf=0 do nn=nharm,nharm+4 ygn=nn*yg if(iwarm.eq.1) then rdu2=2.0d0*(ygn-ygc) u1=anpl+sqrt(rdu2) u2=anpl-sqrt(rdu2) uu2=min(u1*u1,u2*u2) argex=amu*uu2/2.0d0 else rdu2=ygn**2-ygc**2 g1=(ygn+anpl*sqrt(rdu2))/dnl g2=(ygn-anpl*sqrt(rdu2))/dnl gg=min(g1,g2) argex=amu*(gg-1.0d0) u1=(ygn*anpl+sqrt(rdu2))/dnl u2=(ygn*anpl-sqrt(rdu2))/dnl end if if(argex.le.xxcr) nhf=nn end do if(nhf.eq.0) return c lrm=ilarm if(lrm.lt.nhf) lrm=nhf call dispersion(lrm) c akim=ak0*anprim ratiovgr=2.0d0*anpr/derdnm c ratiovgr=2.0d0*anpr/derdnm*vgm alpha=2.0d0*akim*ratiovgr if(alpha.lt.0.0d0) then ierr=94 end if c ithn=1 if(lrm.gt.nharm) ithn=2 if(ieccd.gt.0) call eccd(effjcd) return 999 format(12(1x,e12.5)) end c c c subroutine dispersion(lrm) implicit real*8(a-h,o-z) complex*16 cc0,cc2,cc4,rr complex*16 anpr2a,anpra complex*16 anpr,anpr2,ex,ey,ez,den complex*16 epsl(3,3,lrm),sepsl(3,3),e330 complex*16 e11,e22,e12,e13,e23 c complex*16 e33,e21,e31,e32 complex*16 a13,a31,a23,a32,a33 c common/ygyg/yg common/xgxg/xg common/nplr/anpl,anprf common/mode/sox common/warm/iwarm,ilarm common/nprw/anprr,anpri common/epolar/ex,ey,ez common/amut/amu common/imx/imx errnpr=1.0d0 anpr2a=anprf**2 anpl2=anpl*anpl c if (iwarm.eq.1) then call diel_tens_wr(e330,epsl,lrm) else call diel_tens_fr(e330,epsl,lrm) end if c imxx=abs(imx) c loop to disable convergence if failure detected do do i=1,imxx c do j=1,3 do k=1,3 sepsl(k,j)=dcmplx(0.0d0,0.0d0) do ilrm=1,lrm sepsl(k,j)=sepsl(k,j)+epsl(k,j,ilrm)*anpr2a**(ilrm-1) end do end do end do c anpra=sqrt(anpr2a) c e11=sepsl(1,1) e22=sepsl(2,2) e12=sepsl(1,2) a33=sepsl(3,3) a13=sepsl(1,3) a23=sepsl(2,3) a31=a13 a32=-a23 c e33=e330+anpr2a*a33 e13=anpra*a13 e23=anpra*a23 c e21=-e12 c e31=e13 c e32=-e23 c cc4=(e11-anpl2)*(1.0d0-a33)+(a13+anpl)*(a31+anpl) cc2=-e12*e12*(1.0d0-a33) . -a32*e12*(a13+anpl)+a23*e12*(a31+anpl) . -(a23*a32+e330+(e22-anpl2)*(1.0d0-a33))*(e11-anpl2) . -(a13+anpl)*(a31+anpl)*(e22-anpl2) cc0=e330*((e11-anpl2)*(e22-anpl2)+e12*e12) c rr=cc2*cc2-4.0d0*cc0*cc4 c if(yg.gt.1.0d0) then s=sox if(dimag(rr).LE.0.0d0) s=-s else s=-sox if(dble(rr).le.0.d0.and.dimag(rr).ge.0.d0) s=-s end if c anpr2=(-cc2+s*sqrt(rr))/(2.0d0*cc4) c errnpr=abs(1.0d0-abs(anpr2)/abs(anpr2a)) if(i.gt.1.and.errnpr.lt.1.0d-5) exit anpr2a=anpr2 end do if(i.gt.imxx .and. imxx.gt.1) then if (imx.lt.0) then write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,"// ."': convergence disabled.',e12.5)") xg,yg,sqrt(abs(anpr2)),anpl imxx=1 else write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,"// ."': convergence failed.',e12.5)") xg,yg,sqrt(abs(anpr2)),anpl exit end if else exit end if print*,yg,imx,imxx end do c end loop to disable convergence if(sqrt(dble(anpr2)).lt.0.0d0 .or. anpr2.ne.anpr2 . .or. abs(anpr2).eq.huge(one) .or. abs(anpr2).le.tiny(one)) then write(*,"(' X =',f7.4,' Y =',f7.4,"// . "' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(anpr2)) ierr=99 anpr2=(0.0d0,0.0d0) end if c anpr=sqrt(anpr2) anprr=dble(anpr) anpri=dimag(anpr) c ex=dcmplx(0.0d0,0.0d0) ey=dcmplx(0.0d0,0.0d0) ez=dcmplx(0.0d0,0.0d0) c if (abs(anpl).gt.1.0d-6) then den=e12*e23-(e13+anpr*anpl)*(e22-anpr2-anpl2) ey=-(e12*(e13+anpr*anpl)+(e11-anpl2)*e23)/den ez=(e12*e12+(e22-anpr2-anpl2)*(e11-anpl2))/den ez2=abs(ez)**2 ey2=abs(ey)**2 enx2=1.0d0/(1.0d0+ey2+ez2) ex=dcmplx(sqrt(enx2),0.0d0) ez2=ez2*enx2 ey2=ey2*enx2 ez=ez*ex ey=ey*ex else if(sox.lt.0.0d0) then ez=dcmplx(1.0d0,0.0d0) ez2=abs(ez)**2 else ex2=1.0d0/(1.0d0+abs(-e11/e12)**2) ex=sqrt(ex2) ey=-ex*e11/e12 ey2=abs(ey)**2 ez2=0.0d0 end if end if c return 99 format(20(1x,e12.5)) end c c Fully relativistic case c computation of dielectric tensor elements c up to third order in Larmor radius for hermitian part c subroutine diel_tens_fr(e330,epsl,lrm) implicit real*8(a-h,o-z) complex*16 ui complex*16 e330,epsl(3,3,lrm) complex*16 ca11,ca12,ca22,ca13,ca23,ca33 complex*16 cq0p,cq0m,cq1p,cq1m,cq2p parameter(pi=3.14159265358979d0,rpi=1.77245385090552d0) parameter(ui=(0.0d0,1.0d0)) dimension rr(-lrm:lrm,0:2,0:lrm),ri(lrm,0:2,lrm) c common/xgxg/xg common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/resah/anpl2,dnl c common/cri/cr,ci common/warm/iwarm,ilarm c anpl2=anpl**2 dnl=1.0d0-anpl2 c cmxw=1.0d0+15.0d0/(8.0d0*amu)+105.0d0/(128.0d0*amu**2) cr=-amu*amu/(rpi*cmxw) ci=sqrt(2.0d0*pi*amu)*amu**2/cmxw c do l=1,lrm do j=1,3 do i=1,3 epsl(i,j,l)=dcmplx(0.0d0,0.0d0) end do end do end do c if(iwarm.eq.2) call hermitian(rr,lrm) if(iwarm.eq.4) call hermitian_2(rr,lrm) c call antihermitian(ri,lrm) c do l=1,lrm lm=l-1 fal=-0.25d0**l*fact(2*l)/(fact(l)**2*yg**(2*lm)) ca11=(0.0d0,0.0d0) ca12=(0.0d0,0.0d0) ca13=(0.0d0,0.0d0) ca22=(0.0d0,0.0d0) ca23=(0.0d0,0.0d0) ca33=(0.0d0,0.0d0) do is=0,l k=l-is w=(-1.0d0)**k c asl=w/(fact(is+l)*fact(l-is)) bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) c if(is.gt.0) then cq0p=rr(is,0,l)+rr(-is,0,l)+ui*ri(is,0,l) cq0m=rr(is,0,l)-rr(-is,0,l)+ui*ri(is,0,l) cq1p=rr(is,1,l)+rr(-is,1,l)+ui*ri(is,1,l) cq1m=rr(is,1,l)-rr(-is,1,l)+ui*ri(is,1,l) cq2p=rr(is,2,l)+rr(-is,2,l)+ui*ri(is,2,l) else cq0p=rr(is,0,l) cq0m=rr(is,0,l) cq1p=rr(is,1,l) cq1m=rr(is,1,l) cq2p=rr(is,2,l) end if c ca11=ca11+is**2*asl*cq0p ca12=ca12+is*l*asl*cq0m ca22=ca22+bsl*cq0p ca13=ca13+is*asl*cq1m/yg ca23=ca23+l*asl*cq1p/yg ca33=ca33+asl*cq2p/yg**2 end do epsl(1,1,l) = - xg*ca11*fal epsl(1,2,l) = + ui*xg*ca12*fal epsl(2,2,l) = - xg*ca22*fal epsl(1,3,l) = - xg*ca13*fal epsl(2,3,l) = - ui*xg*ca23*fal epsl(3,3,l) = - xg*ca33*fal end do c cq2p=rr(0,2,0) e330=1.0d0+xg*cq2p c epsl(1,1,1) = 1.d0 + epsl(1,1,1) epsl(2,2,1) = 1.d0 + epsl(2,2,1) c do l=1,lrm epsl(2,1,l) = - epsl(1,2,l) epsl(3,1,l) = epsl(1,3,l) epsl(3,2,l) = - epsl(2,3,l) end do c return end c c c subroutine hermitian(rr,lrm) implicit real*8(a-h,o-z) parameter(tmax=5,npts=500) dimension rr(-lrm:lrm,0:2,0:lrm) real*8 ttv(npts+1),extv(npts+1) c common/ttex/ttv,extv c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/warm/iwarm,ilarm common/cri/cr,ci c do n=-lrm,lrm do k=0,2 do m=0,lrm rr(n,k,m)=0.0d0 end do end do end do c llm=min(3,lrm) c dt=2.0d0*tmax/dble(npts) bth2=2.0d0/amu bth=sqrt(bth2) amu2=amu*amu amu4=amu2*amu2 amu6=amu4*amu2 c do i = 1, npts+1 t = ttv(i) rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) x = t*rxt upl2=bth2*x**2 upl=bth*x gx=1.0d0+t*t/amu exdx=cr*extv(i)*gx/rxt*dt c n1=1 if(iwarm.gt.2) n1=-llm c do n=n1,llm nn=abs(n) gr=anpl*upl+n*yg zm=-amu*(gx-gr) s=amu*(gx+gr) zm2=zm*zm zm3=zm2*zm call calcei3(zm,fe0m) c do m=nn,llm if(n.eq.0.and.m.eq.0) then rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2 else ffe=0.0d0 if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))/amu2 if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s . +s*s*(1.0d0+zm-zm2*fe0m))/amu4 if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0* . (20.0d0-8.0d0*zm+zm2)+s**3*(2.0d0+zm+zm2-zm3*fe0m))/amu6 c rr(n,0,m) = rr(n,0,m) + exdx*ffe rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2 end if c end do end do end do c if(iwarm.gt.2) return c sy1=1.0d0+yg sy2=1.0d0+yg*2.0d0 sy3=1.0d0+yg*3.0d0 c bth4=bth2*bth2 bth6=bth4*bth2 c anpl2=anpl*anpl c rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2) . +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2)) rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2)) rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2)) rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2 . /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/ . sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3)) rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ . 1.5d0*anpl2/sy1**2)) rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* . anpl2/sy1**2)) c if(llm.gt.1) then c rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+ . bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2)) rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2)) rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* . (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4* . (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2 . -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* . (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2)) rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* . (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4* . (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2 . -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* . (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2)) c if(llm.gt.2) then c rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4* . (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2)) rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2)) rr(-1,0,3) = -12.0d0*bth4/sy1* . (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+ . bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2) . /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,3) = -6.0d0*bth6/sy1* . (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2)) c rr(-2,0,3) = -12.0d0*bth4/sy2* . (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+ . bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2) . /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,3) = -6.0d0*bth6/sy2* . (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2)) c rr(-3,0,3) = -12.0d0*bth4/sy3* . (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+ . bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2) . /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4)) rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2)) rr(-3,2,3) = -6.0d0*bth6/sy3* . (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2)) c end if end if c return end subroutine hermitian_2(rr,lrm) implicit real*8(a-h,o-z) parameter(tmax=5,npts=500) dimension rr(-lrm:lrm,0:2,0:lrm) parameter(epsa=0.0d0,epsr=1.0d-4) parameter (lw=5000,liw=lw/4) dimension w(lw),iw(liw) external fhermit c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/warm/iwarm,ilarm common/cri/cr,ci common/nmhermit/n,m,ih c do n=-lrm,lrm do k=0,2 do m=0,lrm rr(n,k,m)=0.0d0 end do end do end do c llm=min(3,lrm) c bth2=2.0d0/amu bth=sqrt(bth2) amu2=amu*amu amu4=amu2*amu2 amu6=amu4*amu2 c n1=1 if(iwarm.gt.10) n1=-llm c do n=n1,llm nn=abs(n) do m=nn,llm if(n.eq.0.and.m.eq.0) then ih=2 c call dqagi(fhermit,bound,2,epsa,epsr,resfh, call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, . epp,neval,ier,liw,lw,last,iw,w) rr(n,2,m) = resfh else do ih=0,2 c call dqagi(fhermit,bound,2,epsa,epsr,resfh, call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, . epp,neval,ier,liw,lw,last,iw,w) rr(n,ih,m) = resfh end do end if end do end do if(iwarm.gt.10) return c sy1=1.0d0+yg sy2=1.0d0+yg*2.0d0 sy3=1.0d0+yg*3.0d0 c bth4=bth2*bth2 bth6=bth4*bth2 c anpl2=anpl*anpl c rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2) . +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2)) rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2)) rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2)) rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2 . /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/ . sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3)) rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ . 1.5d0*anpl2/sy1**2)) rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* . anpl2/sy1**2)) c if(llm.gt.1) then c rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+ . bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2)) rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2)) rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* . (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4* . (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2 . -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* . (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2)) rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* . (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4* . (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2 . -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* . (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2)) c if(llm.gt.2) then c rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4* . (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2)) rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2)) rr(-1,0,3) = -12.0d0*bth4/sy1* . (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+ . bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2) . /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,3) = -6.0d0*bth6/sy1* . (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2)) c rr(-2,0,3) = -12.0d0*bth4/sy2* . (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+ . bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2) . /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,3) = -6.0d0*bth6/sy2* . (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2)) c rr(-3,0,3) = -12.0d0*bth4/sy3* . (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+ . bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2) . /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4)) rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2)) rr(-3,2,3) = -6.0d0*bth6/sy3* . (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2)) c end if end if c return end function fhermit(t) implicit real*8 (a-h,o-z) c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/cri/cr,ci common/nmhermit/n,m,ih bth2=2.0d0/amu bth=sqrt(bth2) amu2=amu*amu amu4=amu2*amu2 amu6=amu4*amu2 rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) x = t*rxt upl2=bth2*x**2 upl=bth*x gx=1.0d0+t*t/amu exdxdt=cr*exp(-t*t)*gx/rxt nn=abs(n) gr=anpl*upl+n*yg zm=-amu*(gx-gr) s=amu*(gx+gr) zm2=zm*zm zm3=zm2*zm call calcei3(zm,fe0m) ffe=0.0d0 uplh=upl**ih if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/amu2 if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s . +s*s*(1.0d0+zm-zm2*fe0m))*uplh/amu4 if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2) . +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/amu6 fhermit= exdxdt*ffe return end c c c subroutine antihermitian(ri,lrm) implicit none integer lmx,nmx,lrm,n,k,m,mm real*8 rpi parameter(rpi=1.77245385090552d0) parameter(lmx=20,nmx=lmx+2) real*8 fsbi(nmx) real*8 ri(lrm,0:2,lrm) real*8 yg,amu,anpl,cmu,anpl2,dnl,cr,ci,ygn,rdu2,rdu,anpr real*8 du,ub,aa,up,um,gp,gm,xp,xm,eep,eem,ee,cm,cim real*8 fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0 real*8 fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m real*8 fact c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/uu/ub,cmu common/resah/anpl2,dnl c common/cri/cr,ci c do n=1,lrm do k=0,2 do m=1,lrm ri(n,k,m)=0.0d0 end do end do end do c dnl=1.0d0-anpl2 cmu=anpl*amu c do n=1,lrm ygn=n*yg rdu2=ygn**2-dnl if(rdu2.gt.0.0d0) then rdu=sqrt(rdu2) du=rdu/dnl ub=anpl*ygn/dnl aa=amu*anpl*du if (abs(aa).gt.5.0d0) then up=ub+du um=ub-du gp=anpl*up+ygn gm=anpl*um+ygn xp=up+1.0d0/cmu xm=um+1.0d0/cmu eem=exp(-amu*(gm-1.0d0)) eep=exp(-amu*(gp-1.0d0)) fi0p0=-1.0d0/cmu fi1p0=-xp/cmu fi2p0=-(1.0d0/cmu**2+xp*xp)/cmu fi0m0=-1.0d0/cmu fi1m0=-xm/cmu fi2m0=-(1.0d0/cmu**2+xm*xm)/cmu do m=1,lrm fi0p1=-2.0d0*m*(fi1p0-ub*fi0p0)/cmu fi0m1=-2.0d0*m*(fi1m0-ub*fi0m0)/cmu fi1p1=-((1.0d0+2.0d0*m)*fi2p0-2.0d0*(m+1.0d0)*ub*fi1p0 . +up*um*fi0p0)/cmu fi1m1=-((1.0d0+2.0d0*m)*fi2m0-2.0d0*(m+1.0d0)*ub*fi1m0 . +up*um*fi0m0)/cmu fi2p1=(2.0d0*(1.0d0+m)*fi1p1-2.0d0*m* . (ub*fi2p0-up*um*fi1p0))/cmu fi2m1=(2.0d0*(1.0d0+m)*fi1m1-2.0d0*m* . (ub*fi2m0-up*um*fi1m0))/cmu if(m.ge.n) then ri(n,0,m)=0.5d0*ci*dnl**m*(fi0p1*eep-fi0m1*eem) ri(n,1,m)=0.5d0*ci*dnl**m*(fi1p1*eep-fi1m1*eem) ri(n,2,m)=0.5d0*ci*dnl**m*(fi2p1*eep-fi2m1*eem) end if fi0p0=fi0p1 fi1p0=fi1p1 fi2p0=fi2p1 fi0m0=fi0m1 fi1m0=fi1m1 fi2m0=fi2m1 end do else ee=exp(-amu*(ygn-1.0d0+anpl*ub)) call ssbi(aa,n,lrm,fsbi) do m=n,lrm cm=rpi*fact(m)*du**(2*m+1) cim=0.5d0*ci*dnl**m mm=m-n+1 fi0m=cm*fsbi(mm) fi1m=-0.5d0*aa*cm*fsbi(mm+1) fi2m=0.5d0*cm*(fsbi(mm+1)+0.5d0*aa*aa*fsbi(mm+2)) ri(n,0,m)=cim*ee*fi0m ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m) ri(n,2,m)=cim*ee*(du*du*fi2m+2.0d0*du*ub*fi1m+ub*ub*fi0m) end do end if end if end do c return end c subroutine ssbi(zz,n,l,fsbi) implicit none integer n,l,nmx,lmx,k,m,mm real*8 eps,c0,c1,sbi,zz real*8 gamm parameter(eps=1.0d-10,lmx=20,nmx=lmx+2) real*8 fsbi(nmx) do m=n,l+2 c0=1.0d0/gamm(dble(m)+1.5d0) sbi=c0 do k=1,50 c1=c0*0.25d0*zz**2/(dble(m+k)+0.5d0)/dble(k) sbi=sbi+c1 if(c1/sbi.lt.eps) exit c0=c1 end do mm=m-n+1 fsbi(mm)=sbi end do return end c c Weakly relativistic dielectric tensor c computation of dielectric tensor elements c Krivenki and Orefice, JPP 30,125 (1983) c subroutine diel_tens_wr(ce330,cepsl,lrm) implicit real*8 (a-b,d-h,o-z) implicit complex*16 (c) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm) parameter(cui=(0.0d0,1.0d0)) c common/xgxg/xg common/ygyg/yg common/nplr/anpl,anprf common/amut/amu c anpl2=anpl*anpl c call fsup(cefp,cefm,lrm) c do l=1,lrm lm=l-1 fcl=0.5d0**l*((1.0d0/yg)**2/amu)**lm*fact(2*l)/fact(l) ca11=(0.d0,0.d0) ca12=(0.d0,0.d0) ca13=(0.d0,0.d0) ca22=(0.d0,0.d0) ca23=(0.d0,0.d0) ca33=(0.d0,0.d0) do is=0,l k=l-is w=(-1.0d0)**k c asl=w/(fact(is+l)*fact(l-is)) bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) c cq0p=amu*cefp(is,0) cq0m=amu*cefm(is,0) cq1p=amu*anpl*(cefp(is,0)-cefp(is,1)) cq1m=amu*anpl*(cefm(is,0)-cefm(is,1)) cq2p=cefp(is,1)+amu*anpl2*(cefp(is,2) . +cefp(is,0)-2.0d0*cefp(is,1)) c ca11=ca11+is**2*asl*cq0p ca12=ca12+is*l*asl*cq0m ca22=ca22+bsl*cq0p ca13=ca13+is*asl*cq1m/yg ca23=ca23+l*asl*cq1p/yg ca33=ca33+asl*cq2p/yg**2 end do cepsl(1,1,l) = - xg*ca11*fcl cepsl(1,2,l) = + cui*xg*ca12*fcl cepsl(2,2,l) = - xg*ca22*fcl cepsl(1,3,l) = - xg*ca13*fcl cepsl(2,3,l) = - cui*xg*ca23*fcl cepsl(3,3,l) = - xg*ca33*fcl end do c cq2p=cefp(0,1)+amu*anpl2*(cefp(0,2)+cefp(0,0)-2.0d0*cefp(0,1)) ce330=1.0d0-xg*amu*cq2p c cepsl(1,1,1) = 1.d0 + cepsl(1,1,1) cepsl(2,2,1) = 1.d0 + cepsl(2,2,1) c do l=1,lrm cepsl(2,1,l) = - cepsl(1,2,l) cepsl(3,1,l) = cepsl(1,3,l) cepsl(3,2,l) = - cepsl(2,3,l) end do c return end c c c subroutine fsup(cefp,cefm,lrm) implicit real*8 (a-b,d-h,o-z) implicit complex*16 (c) parameter(apsicr=0.7d0) parameter(cui=(0.0d0,1.0d0)) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2) c common/ygyg/yg common/nplr/anpl,anprf common/amut/amu c psi=sqrt(0.5d0*amu)*anpl apsi=abs(psi) c do is=0,lrm alpha=anpl*anpl/2.0d0+is*yg-1.0d0 phi2=amu*alpha phim=sqrt(abs(phi2)) if (alpha.ge.0) then xp=psi-phim yp=0.0d0 xm=-psi-phim ym=0.0d0 x0=-phim y0=0.0d0 else xp=psi yp=phim xm=-psi ym=phim x0=0.0d0 y0=phim end if call zetac (xp,yp,zrp,zip,iflag) call zetac (xm,ym,zrm,zim,iflag) c czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) cf12=(0.0d0,0.0d0) if (alpha.ge.0) then if (alpha.ne.0) cf12=-(czp+czm)/(2.0d0*phim) else cf12=-cui*(czp+czm)/(2.0d0*phim) end if c if(apsi.gt.apsicr) then cf32=-(czp-czm)/(2.0d0*psi) else cphi=phim if(alpha.lt.0) cphi=-cui*phim call zetac (x0,y0,zr0,zi0,iflag) cz0=dcmplx(zr0,zi0) cdz0=2.0d0*(1.0d0-cphi*cz0) cf32=cdz0 end if c cf0=cf12 cf1=cf32 cefp(is,0)=cf32 cefm(is,0)=cf32 do l=1,is+2 iq=l-1 if(apsi.gt.apsicr) then cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 else cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) end if ir=l-is if(ir.ge.0) then cefp(is,ir)=cf2 cefm(is,ir)=cf2 end if cf0=cf1 cf1=cf2 end do c if(is.ne.0) then c alpha=anpl*anpl/2.0d0-is*yg-1.0d0 phi2=amu*alpha phim=sqrt(abs(phi2)) if (alpha.ge.0.0d0) then xp=psi-phim yp=0.0d0 xm=-psi-phim ym=0.0d0 x0=-phim y0=0.0d0 else xp=psi yp=phim xm=-psi ym=phim x0=0.0d0 y0=phim end if call zetac (xp,yp,zrp,zip,iflag) call zetac (xm,ym,zrm,zim,iflag) c czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) c cf12=(0.0d0,0.0d0) if (alpha.ge.0) then if (alpha.ne.0.0d0) cf12=-(czp+czm)/(2.0d0*phim) else cf12=-cui*(czp+czm)/(2.0d0*phim) end if if(apsi.gt.apsicr) then cf32=-(czp-czm)/(2.0d0*psi) else cphi=phim if(alpha.lt.0) cphi=-cui*phim call zetac (x0,y0,zr0,zi0,iflag) cz0=dcmplx(zr0,zi0) cdz0=2.0d0*(1.0d0-cphi*cz0) cf32=cdz0 end if c cf0=cf12 cf1=cf32 do l=1,is+2 iq=l-1 if(apsi.gt.apsicr) then cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 else cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) end if ir=l-is if(ir.ge.0) then cefp(is,ir)=cefp(is,ir)+cf2 cefm(is,ir)=cefm(is,ir)-cf2 end if cf0=cf1 cf1=cf2 end do c end if c end do c return end subroutine eccd(effjcd) use green_func_p implicit none real*8 anum,denom,resp,resj,effjcd,coullog,dens,tekev real*8 anucc,canucc,ddens real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc real*8 fjch,fjncl,fjch0,fconic real*8 cst,cst2 integer ieccd,nharm,nhf,nhn external fjch,fjncl,fjch0 parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) parameter(vcsi=2.99792458d+8) parameter(qe=4.8032068d-10,me=mesi*1.d3,vc=vcsi*1.d2) parameter(canucc=4.0d13*pi*qe**4/(me**2*vc**3)) parameter(ceff=qesi/(mesi*vcsi)) common/nharm/nharm,nhf common/nhn/nhn common/ieccd/ieccd common/tete/tekev common/dddens/dens,ddens common/zz/Zeff common/btot/btot common/bmxmn/bmax,bmin common/fc/fc common/ncl/rbx common/cohen/rbn,alams,pa,fp0s common/cst/cst,cst2 anum=0.0d0 denom=0.0d0 effjcd=0.0d0 coullog=24.0d0-log(1.0d4*sqrt(0.1d0*dens)/tekev) anucc=canucc*dens*coullog c nhf=nharm select case(ieccd) case(1) c cohen model c rbn=B/B_min c rbx=B/B_max c cst2=1.0d0-B/B_max c alams=sqrt(1-B_min/B_max) c Zeff < 31 !!! c fp0s= P_a (alams) rbn=btot/bmin rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 cst=sqrt(cst2) alams=sqrt(1.0d0-bmin/bmax) pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0 fp0s=fconic(alams,pa,0) do nhn=nharm,nhf call curr_int(fjch,resj,resp) anum=anum+resj denom=denom+resp end do case(2:9) cst=0.0d0 cst2=0.0d0 do nhn=nharm,nhf call curr_int(fjch0,resj,resp) anum=anum+resj denom=denom+resp end do case(10:11) c neoclassical model: c ft=1-fc trapped particle fraction c rzfc=(1+Zeff)/fc rbn=btot/bmin rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 cst=sqrt(cst2) call SpitzFuncCoeff(Tekev,Zeff,fc) do nhn=nharm,nhf call curr_int(fjncl,resj,resp) anum=anum+resj denom=denom+resp end do nhn=nhn-1 CASE DEFAULT write(*,*) 'ieccd undefined' end select c c effjpl = / /(B_min/) [A m /W] c if(denom.gt.0.0d0) effjcd=-ceff*anum/(anucc*denom) return end c c c subroutine curr_int(fcur,resj,resp) implicit real*8(a-h,o-z) parameter(epsa=0.0d0,epsr=1.0d-2) parameter(xxcr=16.0d0) parameter (lw=5000,liw=lw/4) dimension w(lw),iw(liw) external fcur,fpp common/nhn/nhn common/ygyg/yg common/nplr/anpl,anpr common/amut/amu common/gg/uplp,uplm,ygn common/ierr/ierr common/iokh/iokhawa common/cst/cst,cst2 c EC power and current densities anpl2=anpl*anpl dnl=1.0d0-anpl2 ygn=nhn*yg ygn2=ygn*ygn resj=0.0d0 resj1=0.0d0 resj2=0.0d0 resp=0.0d0 c rdu2=anpl2+ygn2-1.0d0 uplp=0.0d0 uplm=0.0d0 upltp=0.0d0 upltm=0.0d0 c if (rdu2.ge.0.0d0) then rdu=sqrt(rdu2) uplp=(anpl*ygn+rdu)/dnl uplm=(anpl*ygn-rdu)/dnl c uu1=uplm uu2=uplp xx1=amu*(anpl*uu1+ygn-1.0d0) xx2=amu*(anpl*uu2+ygn-1.0d0) c if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl duu=abs(uu1-uu2) c if(duu.gt.1.d-6) then call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ier, . liw,lw,last,iw,w) if (ier.gt.0) ierr=90 end if c rdu2t=cst2*anpl2+ygn2-1.0d0 c if (rdu2t.lt.0.0d0.or.cst2.eq.0.0d0) then c c resonance curve does not cross the trapping region c iokhawa=0 if(duu.gt.1.d-4) then call dqags(fcur,uu1,uu2,epsa,epsr, . resj,ej,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) ierr=91 end if else c c resonance curve crosses the trapping region c iokhawa=1 rdut=sqrt(rdu2t) upltm=(cst2*anpl*ygn-cst*rdut)/(1.0d0-cst2*anpl2) upltp=(cst2*anpl*ygn+cst*rdut)/(1.0d0-cst2*anpl2) c uu1=uplm uu2=upltm xx1=amu*(anpl*uu1+ygn-1.0d0) xx2=amu*(anpl*uu2+ygn-1.0d0) if(xx1.lt.xxcr.or.xx2.lt.xxcr) then if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl duu=abs(uu1-uu2) if(duu.gt.1.d-6) then call dqags(fcur,uu1,uu2,epsa,epsr, . resj1,ej1,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then if (abs(resj1).lt.1.0d-10) then resj1=0.0d0 else ierr=92 end if end if end if end if c uu1=upltp uu2=uplp xx1=amu*(anpl*uu1+ygn-1.0d0) xx2=amu*(anpl*uu2+ygn-1.0d0) if(xx1.lt.xxcr.or.xx2.lt.xxcr) then if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl duu=abs(uu1-uu2) if(duu.gt.1.d-6) then call dqags(fcur,uu1,uu2,epsa,epsr, . resj2,ej2,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then if(ier.ne.2) ierr=93 end if end if end if resj=resj1+resj2 end if end if return end c c computation of integral for power density, integrand function fpp c c ith=0 : polarization term = const c ith=1 : polarization term Larmor radius expansion to lowest order c ith=2 : full polarization term (J Bessel) c function fpp(upl) implicit real*8 (a-h,o-z) complex*16 ex,ey,ez,ui,emxy,epxy parameter(ui=(0.0d0,1.0d0)) common/ygyg/yg common/nplr/anpl,anpr common/amut/amu common/gg/uplp,uplm,ygn common/epolar/ex,ey,ez common/nprw/anprre,anprim common/ithn/ithn common/nhn/nhn upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm) gam=anpl*upl+ygn ee=exp(-amu*(gam-1)) thn2=1.0d0 thn2u=upr2*thn2 if(ithn.gt.0) then emxy=ex-ui*ey epxy=ex+ui*ey if(upr2.gt.0.0d0) then upr=sqrt(upr2) bb=anprre*upr/yg if(ithn.eq.1) then c Larmor radius expansion polarization term at lowest order cth=1.0d0 if(nhn.gt.1) cth=(0.5d0*bb)**(nhn-1)*nhn/fact(nhn) thn2=(0.5d0*cth*abs(emxy+ez*anprre*upl/ygn))**2 thn2u=upr2*thn2 else c Full polarization term nm=nhn-1 np=nhn+1 ajbnm=dbesjn(nm, bb) ajbnp=dbesjn(np, bb) ajbn=dbesjn(nhn, bb) thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0d0))**2 end if end if end if fpp=ee*thn2u return end c computation of integral for current density c fjch integrand for Cohen model with trapping c fjch0 integrand for Cohen model without trapping c fjncl integrand for momentum conserv. model K(u) from Maruschenko c gg=F(u)/u with F(u) as in Cohen paper function fjch(upl) implicit real*8 (a-h,o-z) common/gg/uplp,uplm,ygn common/nplr/anpl,anpr common/zz/Zeff common/cohen/rb,alams,pa,fp0s upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm) gam=anpl*upl+ygn u2=gam*gam-1.0d0 u=sqrt(u2) z5=Zeff+5.0d0 xi=1.0d0/z5**2 xib=1.0d0-xi xibi=1.0d0/xib fu2b=1.0d0+xib*u2 fu2=1.0d0+xi*u2 gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2) gg=u*gu/z5 dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 alam=sqrt(1.0d0-upr2/u2/rb) fp0=fconic(alam,pa,0) dfp0=-(pa*pa/2.0d0+0.125d0) if (alam.lt.1.0d0) then dfp0=-fconic(alam,pa,1)/sqrt(1.0d0-alam**2) end if fh=alam*(1.0d0-alams*fp0/(alam*fp0s)) dfhl=1.0d0-alams*dfp0/fp0s eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam) if(upl.lt.0.0d0) eta=-eta fjch=eta*fpp(upl) return end function fjch0(upl) implicit real*8 (a-h,o-z) common/gg/uplp,uplm,ygn common/nplr/anpl,anpr common/zz/Zeff gam=anpl*upl+ygn u2=gam*gam-1.0d0 u=sqrt(u2) z5=Zeff+5.0d0 xi=1.0d0/z5**2 xib=1.0d0-xi xibi=1.0d0/xib fu2b=1.0d0+xib*u2 fu2=1.0d0+xi*u2 gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2) gg=u*gu/z5 dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 eta=anpl*gg+gam*upl*dgg/u fjch0=eta*fpp(upl) return end function fjncl(upl) use green_func_p implicit real*8 (a-h,o-z) common/gg/uplp,uplm,ygn common/nplr/anpl,anpr common/fc/fc common/ncl/hb common/psinv/psinv common/amut/amu common/tete/tekev common/zz/Zeff gam=anpl*upl+ygn u2=gam*gam-1.0d0 u=sqrt(u2) upr2=u2-upl*upl bth=sqrt(2.0d0/amu) uth=u/bth call GenSpitzFunc(Tekev,Zeff,fc,uth,u,gam,fk,dfk) fk=fk*(4.0d0/amu**2) dfk=dfk*(2.0d0/amu)*bth alam=upr2/u2/hb psi=psinv call vlambda(alam,psi,fu,dfu) eta=gam*fu*dfk/u-2.0d0*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb if(upl.lt.0) eta=-eta fjncl=eta*fpp(upl) return end subroutine vlambda(alam,psi,fv,dfv) implicit real*8 (a-h,o-z) parameter (nnintp=101,nlam=41) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ . njest*nnintp+nlest+54) parameter(kwrk=nnintp+nlam+njest+nlest+3) parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) external dierckx_fpbisp dimension xxs(1),yys(1),ffs(1) dimension ch01(lw01),ch((njest-4)*(nlest-4)) dimension tjp(njest),tlm(nlest) dimension iwrk(kwrk),wrk(lwrk) common/coffvl/ch common/coffdvl/ch01 common/coffvlt/tjp,tlm common/coffvln/njpt,nlmt njp=njpt nlm=nlmt xxs(1)=sqrt(psi) yys(1)=alam call dierckx_fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs, . wrk(1),wrk(5),iwrk(1),iwrk(2)) fv=ffs(1) iwp=1+(njp-4)*(nlm-5) iwl=iwp+4 call dierckx_fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2, . xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2)) dfv=ffs(1) return end subroutine projxyzt(iproj,nfile) use beamdata, only : ywrk,nrayr,nrayth implicit real*8 (a-h,o-z) c common/psinv11/psinv11 common/istep/istep common/ss/st c rtimn=1.d+30 rtimx=-1.d-30 c jd=1 if(iproj.eq.0) jd=nrayr-1 do j=1,nrayr,jd kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk c dx=ywrk(1,j,k)-ywrk(1,1,1) dy=ywrk(2,j,k)-ywrk(2,1,1) dz=ywrk(3,j,k)-ywrk(3,1,1) c dirx=ywrk(4,j,k) diry=ywrk(5,j,k) dirz=ywrk(6,j,k) dir=sqrt(dirx*dirx+diry*diry+dirz*dirz) c if(j.eq.1.and.k.eq.1) then csth1=dirz/dir snth1=sqrt(1.0d0-csth1**2) csps1=1.0d0 snps1=0.0d0 if(snth1.gt.0.0d0) then csps1=diry/(dir*snth1) snps1=dirx/(dir*snth1) end if end if xti= dx*csps1-dy*snps1 yti=(dx*snps1+dy*csps1)*csth1-dz*snth1 zti=(dx*snps1+dy*csps1)*snth1+dz*csth1 rti=sqrt(xti**2+yti**2) c dirxt= (dirx*csps1-diry*snps1)/dir diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir c if(k.eq.1) then xti1=xti yti1=yti zti1=zti rti1=rti end if if(.not.(iproj.eq.0.and.j.eq.1)) . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 c if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti c end do c c if(.not.(iproj.eq.0.and.j.eq.1)) c . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 if(iproj.eq.1) write(nfile,*) ' ' end do c write(nfile,*) ' ' c write(612,99) istep,st,psinv11,rtimn,rtimx return 99 format(i5,12(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) end c c c subroutine pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout) use const_and_precisions, only : pi use beamdata, only : psjki,iiv,ppabs,ccci,pdjki, . nrayr,nrayth,nstep implicit real*8(a-h,o-z) parameter(rtbc=1.0d0) dimension rhopin(nrho),dpdvout(nrho),ajcdout(nrho) real*8, dimension(:), allocatable :: xxi,ypt,yamp real*8, dimension(:), allocatable :: rtab,rhotv,rtab1,ajphiv, . dpdv,dvolt,darea,ratjav,ratjbv,ratjplv,ajplv,ajcdav, . ajcdbv,pins,currins,fi parameter(llmx=21) dimension isev(llmx) c common/istep/istep common/dsds/dst common/ipec/ipec,nnd common/ieccd/ieccd common/index_rt/index_rt c common/cent/btrcen,rcen common/angles/alpha0,beta0 common/iieq/iequil common/parban/b0,rr0m,zr0m,rpam common/taumnx/taumn,taumx,pabstot,currtot common/jmxmn/rhot1,rhot2,aj1,aj2 c common/polcof/psipol,chipol c allocate(xxi(nstep),ypt(nstep),yamp(nstep)) stf=istep*dst if (ipec.ge.0) then nd=nnd else ! ipec=-1 uses rho_pol input grid nd=nrho end if allocate(rtab(nd),rhotv(nd),rtab1(0:nd),ajphiv(nd), . dpdv(nd),dvolt(nd),darea(nd),ratjav(nd), . ratjbv(nd),ratjplv(nd),ajplv(nd),ajcdav(nd), . ajcdbv(nd),pins(nd),currins(nd),fi(nd)) if(pabstot.gt.0.0d0) then do ll=1,llmx isev(ll)=0 end do intp=0 voli0=0.0d0 areai0=0.0d0 rtab1(0)=0.0d0 do it=1,nd if (ipec.ge.0) then drt=1.0d0/dble(nd-1) rt=dble(it-1)*drt if(it.lt.nd) then rt1=rt+drt/2.0d0 else rt1=rt end if else c radial coordinate of i-(i+1) interval left point rt=rhopin(it) if(it.lt.nd) then c radial coordinate of i-(i+1) interval mid point rt1=(rt+rhopin(it+1))/2.0d0 else rt1=rt end if end if rtab(it)=rt rtab1(it)=rt1 dpdv(it)=0.0d0 ajphiv(it)=0.0d0 if (ipec.eq.0) then psit=rt psit1=rt1 else psit=rt**2 psit1=rt1**2 end if if(iequil.eq.2) then rhotv(it)=frhotor(psit) call valpsispl(sqrt(psit1),voli1,dervoli,areai1,rrii, . rbavi,bmxi,bmni,fci,intp) dvolt(it)=abs(voli1-voli0) darea(it)=abs(areai1-areai0) voli0=voli1 areai0=areai1 call ratioj(sqrt(psit),ratjai,ratjbi,ratjpli) ratjav(it)=ratjai ratjbv(it)=ratjbi ratjplv(it)=ratjpli else rhotv(it)=sqrt(psit) area1=pi*psit1*rpam**2 voli1=2.0d0*pi*rr0m*area1 dvolt(it)=abs(voli1-voli0) darea(it)=abs(areai1-areai0) voli0=voli1 areai0=areai1 end if end do kkk=1 do j=1,nrayr if(j.gt.1) kkk=nrayth do k=1,kkk ise0=0 ii=iiv(j,k) if (ii.lt.nstep) then if(psjki(j,k,ii+1).ne.0.0d0) ii=ii+1 end if idecr=-1 is=1 do i=1,ii if(ipec.eq.0) then xxi(i)=abs(psjki(j,k,i)) else xxi(i)=sqrt(abs(psjki(j,k,i))) end if if(psjki(j,k,i).ge.0.and.psjki(j,k,i).le.rtbc) then ypt(i)=ppabs(j,k,i) yamp(i)=ccci(j,k,i) else ypt(i)=0.0d0 yamp(i)=0.0d0 end if if(ise0.eq.0) then if(xxi(i).lt.rtbc) then ise0=i isev(is)=max(i-1,1) is=is+1 end if else if (idecr.eq.-1) then if(xxi(i).gt.xxi(i-1)) then isev(is)=i-1 is=is+1 idecr=1 end if else if(xxi(i).gt.rtbc) exit if(xxi(i).lt.xxi(i-1)) then ! isev(is)=i !!!!!!!!!! it should be isev(is)=i-1 isev(is)=i-1 is=is+1 idecr=-1 end if end if end if end do isev(is)=max(i-1,1) c ppa1=0.0d0 cci1=0.0d0 do iis=1,is-1 iis1=iis+1 iise0=isev(iis) iise=isev(iis1) if (mod(iis,2).ne.0) then idecr=-1 ind1=nd ind2=2 iind=-1 else idecr=1 ind1=1 ind2=nd iind=1 end if do ind=ind1,ind2,iind !!!!!!!!!! is it safe? ! iind=ind !!!!!!!!!! iind reused in the loop! indi=ind !!!!!!!!!! iind reused in the loop! if (idecr.eq.-1) indi=ind-1 rt1=rtab1(indi) call locatex(xxi,iise,iise0,iise,rt1,itb1) if(itb1.ge.iise0.and.itb1.lt.iise) then call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1), . ypt(itb1+1),rt1,ppa2) call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1), . yamp(itb1+1),rt1,cci2) dppa=ppa2-ppa1 dpdv(ind)=dpdv(ind)+dppa didst=(cci2-cci1) ajphiv(ind)=ajphiv(ind)+didst ppa1=ppa2 cci1=cci2 end if end do end do end do end do rhotpav=0.0d0 rhot2pav=0.0d0 rhotjav=0.0d0 rhotjava=0.0d0 rhot2java=0.0d0 c here dpdv is still dP (not divided yet by dV) spds=sum(dpdv(1:nd)) if (spds.gt.0.0d0) then rhotpav=sum(rhotv(1:nd)*dpdv(1:nd))/spds rhot2pav=sum(rhotv(1:nd)*rhotv(1:nd)*dpdv(1:nd))/spds end if c here ajphiv is still dI (not divided yet by dA) if (ieccd.ne.0) then sccs=sum(ajphiv(1:nd)) if (sccs.gt.0.0d0) then rhotjav=sum(rhotv(1:nd)*ajphiv(1:nd))/sccs end if sccsa=sum(abs(ajphiv(1:nd))) if (sccsa.gt.0.0d0) then rhotjava=sum(rhotv(1:nd)*abs(ajphiv(1:nd)))/sccsa rhot2java=sum(rhotv(1:nd)*rhotv(1:nd)*abs(ajphiv(1:nd)))/sccsa end if end if c factor sqrt(8)=2 sqrt(2) to match with full width c of gaussian profile drhot2pav=rhot2pav-rhotpav**2 drhotpav=sqrt(8.d0*drhot2pav) drhot2java=rhot2java-rhotjava**2 drhotjava=sqrt(8.d0*drhot2java) spds=0.0d0 sccs=0.0d0 do i=1,nd spds=spds+dpdv(i) sccs=sccs+ajphiv(i) pins(i)=spds currins(i)=sccs dpdv(i)=dpdv(i)/dvolt(i) ajphiv(i)=ajphiv(i)/darea(i) end do facpds=1.0d0 facjs=1.0d0 if(spds.gt.0.0d0) facpds=pabs/spds if(sccs.ne.0.0d0) facjs=currt/sccs do i=1,nd dpdv(i)=facpds*dpdv(i) ajphiv(i)=facjs*ajphiv(i) ajcdav(i)=ajphiv(i)*ratjav(i) ajcdbv(i)=ajphiv(i)*ratjbv(i) ajplv(i)=ajphiv(i)*ratjplv(i) end do rhpp=frhopol(rhotpav) rhpj=frhopol(rhotjava) dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, . drhotp,drhopp) if(ieccd.ne.0) then call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi, . drhotjfi,drhopfi) xps=rhopfi c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) else rhotjfi=0.0d0 rhopfi=0.0d0 ajmxfi=0.0d0 ajphip=0.0d0 drhotjfi=0.0d0 drhopfi=0.0d0 xps=rhopp ratjamx=1.0d0 ratjbmx=1.0d0 ratjplmx=1.0d0 end if iif1=iiv(1,1) istmx=1 do i=2,iif1 if(psjki(1,1,i).ge.0.0d0) then if(pdjki(1,1,i).gt.pdjki(1,1,i-1)) istmx=i end if end do stmx=istmx*dst pins_02=0.0d0 pins_05=0.0d0 pins_085=0.0d0 xrhot=0.2d0 call vlocate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_02) xrhot=0.5d0 call vlocate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_05) xrhot=0.85d0 call vlocate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_085) else ajmxfi=0.0d0 ajphip=0.0d0 dpdvp=0.0d0 dpdvmx=0.0d0 rhotjfi=1.0d0 rhotjav=1.0d0 rhotjava=1.0d0 rhotp=1.0d0 rhotpav=1.0d0 drhotjfi=0.0d0 drhotjava=0.0d0 drhotp=0.0d0 drhotpav=0.0d0 ratjamx=1.0d0 ratjbmx=1.0d0 taumn=0 taumx=0 stmx=stf pins_02=0.0d0 pins_05=0.0d0 pins_085=0.0d0 c end of pabstot > 0 end if c dPdV [MW/m^3], Jcd [MA/m^2] if(ieccd.eq.0) currt=0.0d0 currtka=currt*1.0d3 write(*,*)' ' write(*,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' write(*,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp, . rhotjfi,rhotjava,rhotp,rhotpav, . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp write(607,99) currtka,pabstot,ajphip,dpdvp, . rhotjfi,rhotjava,rhotp,rhotpav, . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp do i=1,nd if (ipec.eq.0) then psin=rtab(i) rhop=sqrt(rtab(i)) else psin=rtab(i)**2 rhop=rtab(i) end if pinsr=0.0d0 if(pabstot.gt.0) pinsr=pins(i)/pabstot write(648,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), . currins(i),pins(i),pinsr,real(index_rt) end do c fill output arguments only if ipec=-1 if (ipec.eq.-1) then do i=1,nd dpdvout(i)=dpdv(i) ajcdout(i)=ajcdbv(i) end do end if deallocate(xxi,ypt,yamp) deallocate(rtab,rhotv,rtab1,ajphiv,dpdv,dvolt,darea,ratjav,ratjbv, . ratjplv,ajplv,ajcdav,ajcdbv,pins,currins,fi) return 49 format(i5,20(1x,e12.5)) 99 format(30(1x,e12.5)) end c c c subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) implicit real*8(a-h,o-z) parameter(emn1=0.367879441171442d0) dimension xx(nd),yy(nd) common/jmxmn/rhotp,rhotm,ypkp,ypkm common/ipec/ipec,nnd common/iieq/iequil c call vmaxmini(yy,nd,ymn,ymx,imn,imx) ypk=0.0d0 xmx=xx(imx) xmn=xx(imn) if (abs(ymx).gt.abs(ymn)) then ipk=imx ypkp=ymx xpkp=xmx if(abs(ymn/ymx).lt.1d-2) ymn=0.0d0 ypkm=ymn xpkm=xmn else ipk=imn ypkp=ymn xpkp=xmn if(abs(ymx/ymn).lt.1d-2) ymx=0.0d0 ypkm=ymx xpkm=xmx end if if(xpkp.gt.0.0d0) then xpk=xpkp ypk=ypkp yye=ypk*emn1 call locatex(yy,nd,1,ipk,yye,ie1) if(ie1.gt.0.and.ie1.lt.nd) then call intlin(yy(ie1),xx(ie1),yy(ie1+1),xx(ie1+1),yye,rte1) else rte1=0.0d0 end if call locatex(yy,nd,ipk,nd,yye,ie2) if(ie2.gt.0.and.ie2.lt.nd) then call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2) else rte2=0.0d0 end if else ipk=2 xpk=xx(2) ypk=yy(2) rte1=0.0d0 yye=ypk*emn1 call vlocate(yy,nd,yye,ie2) call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2) end if c ipk=1 if(ymx.ne.0.and.ymn.ne.0) ipk=2 c drhop=0.0d0 drhot=0.0d0 psie1=0.0d0 psie2=1.0d0 rhopmx=1.0d0 rhopmn=0.0d0 if (ie1.gt.0.or.ie2.gt.0) then if(ipec.eq.0) then rhopmx=sqrt(xpkp) rhopmn=sqrt(xpkm) psie2=rte2 psie1=rte1 drhop=sqrt(rte2)-sqrt(rte1) else rhopmx=xpkp rhopmn=xpkm drhop=rte2-rte1 psie2=rte2**2 psie1=rte1**2 end if end if c if(iequil.eq.2) then rhotmx=frhotor(rhopmx**2) rhotmn=frhotor(rhopmn**2) rhote2=frhotor(psie2) rhote1=frhotor(psie1) drhot=rhote2-rhote1 else rhotmx=rhopmx rhotmn=rhopmn drhot=drhop end if c if(ipk.eq.2) then drhop=-drhop drhot=-drhot end if ypk=ypkp rhotp=rhotmx rhotm=rhotmn c return end c subroutine polarcold(exf,eyif,ezf,elf,etf) implicit real*8(a-h,o-z) common/nplr/anpl,anpr common/xgxg/xg common/ygyg/yg common/mode/sox c c dcold dispersion c dielectric tensor (transposed) c c exf=0.0d0 c eyif=0.0d0 c ezf=0.0d0 c if(xg.le.0) return c anpl2=anpl*anpl anpr2=anpr*anpr an2=anpl2+anpr2 yg2=yg**2 dy2=1.0d0-yg2 aa=1.0d0-xg-yg2 e3=1.0d0-xg c if(xg.gt.0.0d0) then if (anpl.ne.0.0d0) then qq=xg*yg/(an2*dy2-aa) p=(anpr2-e3)/(anpl*anpr) ezf=1.0d0/sqrt(1.0d0+p*p*(1.0d0+qq*qq)) exf=p*ezf eyif=qq*exf else if(sox.lt.0.d0) then ezf=1 exf=0 eyif=0 else ezf=0 qq=-aa/(xg*yg) exf=1.0d0/sqrt(1.0d0+qq*qq) eyif=qq*exf end if end if elf=(anpl*ezf+anpr*exf)/sqrt(an2) etf=sqrt(1.0d0-elf*elf) else if(sox.lt.0.0d0) then ezf=1 exf=0 eyif=0 else ezf=0 exf=0.0d0 eyif=1.0d0 end if elf=0 etf=1 end if c return end subroutine pol_limit(qq,uu,vv) use const_and_precisions, only : pi implicit none integer*4 ipolc real*8 bv(3),anv(3) real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2 real*8 deno,denx,anpl2,dnl,del0 real*8 uuom,vvom,qqom,uuxm,vvxm,qqxm,ellom,ellxm,qq,uu,vv real*8 aaom,bbom,llmom,aaxm,bbxm,llmxm,psiom,psixm,chiom,chixm real*8 beta0,alpha0,gam real*8 sox,psipol,chipol complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt parameter(ui=(0.0d0,1.0d0)) c common/anv/anv common/nplr/anpl,anpr common/ygyg/yg common/bb/bv common/angles/alpha0,beta0 common/mode/sox common/polcof/psipol,chipol common/ipol/ipolc common/evt/ext,eyt c anx=anv(1) any=anv(2) anz=anv(3) anpl2=anpl*anpl dnl=1.0d0-anpl2 del0=sqrt(dnl**2+4.0d0*anpl2/yg**2) ffo=0.5d0*yg*(dnl+del0) ffx=0.5d0*yg*(dnl-del0) an2=anx*anx+any*any+anz*anz an=sqrt(an2) anxy=sqrt(anx*anx+any*any) sngam=(anz*anpl-an2*bv(3))/(an*anxy*anpr) csgam=-(any*bv(1)-anx*bv(2))/anxy/anpr csg2=csgam**2 sng2=sngam**2 ffo2=ffo*ffo ffx2=ffx*ffx deno=ffo2+anpl2 denx=ffx2+anpl2 xe2om=(ffo2*csg2+anpl2*sng2)/deno ye2om=(ffo2*sng2+anpl2*csg2)/deno xe2xm=(ffx2*csg2+anpl2*sng2)/denx ye2xm=(ffx2*sng2+anpl2*csg2)/denx c exom=(ffo*csgam-ui*anpl*sngam)/sqrt(deno) eyom=(-ffo*sngam-ui*anpl*csgam)/sqrt(deno) qqom=abs(exom)**2-abs(eyom)**2 uuom=2.0d0*dble(exom*dconjg(eyom)) vvom=2.0d0*dimag(exom*dconjg(eyom)) llmom=sqrt(qqom**2+uuom**2) aaom=sqrt((1+llmom)/2.0d0) bbom=sqrt((1-llmom)/2.0d0) ellom=bbom/aaom psiom=0.5d0*atan2(uuom,qqom)*180.d0/pi chiom=0.5d0*asin(vvom)*180.d0/pi c exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx) eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx) qqxm=abs(exxm)**2-abs(eyxm)**2 uuxm=2.0d0*dble(exxm*dconjg(eyxm)) vvxm=2.0d0*dimag(exxm*dconjg(eyxm)) llmxm=sqrt(qqxm**2+uuxm**2) aaxm=sqrt((1+llmxm)/2.0d0) bbxm=sqrt((1-llmxm)/2.0d0) ellxm=bbxm/aaxm psixm=0.5d0*atan2(uuxm,qqxm)*180.d0/pi chixm=0.5d0*asin(vvxm)*180.d0/pi c if (sox.lt.0.0d0) then psipol=psiom chipol=chiom ext=exom eyt=eyom qq=qqom vv=vvom uu=uuom else psipol=psixm chipol=chixm ext=exxm eyt=eyxm qq=qqxm vv=vvxm uu=uuxm endif gam=atan(sngam/csgam)*180.d0/pi return 111 format(20(1x,e12.5)) end subroutine wall_refl(xvrfl,anvrfl,qqtr,uutr,vvtr,irfl) use const_and_precisions, only : pi implicit none integer*4 ivac,irfl real*8 anv(3),xv(3),xvrfl(3) real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3),xvout(3) real*8 uutr,vvtr,qqtr,qq,uu,vv real*8 psipol,chipol,psitr,chitr complex*16 ui,extr,eytr,eztr,ext,eyt complex*16 evin(3),evrfl(3) parameter(ui=(0.0d0,1.0d0)) common/xv/xv common/anv/anv common/polcof/psipol,chipol common/evt/ext,eyt common/wrefl/walln anv=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) c computation of reflection coordinates and normal to the wall irfl=1 ivac=1 call vacuum_rt(xv,xvout,ivac) if(ivac.lt.0) then irfl=0 xvrfl=xvout xv=xvout anvrfl=anv return end if c rotation matrix from local to lab frame vv1(1)=anv(2) vv1(2)=-anv(1) vv1(3)=0.0d0 vv2(1)=anv(1)*anv(3) vv2(2)=anv(2)*anv(3) vv2(3)=-anv(1)*anv(1)-anv(2)*anv(2) vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2) vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2) vv3=anv evin=ext*vv1+eyt*vv2 c wave vector and electric field after reflection in lab frame anvrfl=anv-2.0d0* . (anv(1)*walln(1)+anv(2)*walln(2)+anv(3)*walln(3))*walln evrfl=-evin+2.0d0* . (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln vv1(1)=anvrfl(2) vv1(2)=-anvrfl(1) vv1(3)=0.0d0 vv2(1)=anvrfl(1)*anvrfl(3) vv2(2)=anvrfl(2)*anvrfl(3) vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2) vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2) vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2) vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2) extr=dot_product(vv1,evrfl) eytr=dot_product(vv2,evrfl) eztr=dot_product(vv3,evrfl) qqtr=abs(extr)**2-abs(eytr)**2 uutr=2.0d0*dble(extr*dconjg(eytr)) vvtr=2.0d0*dimag(extr*dconjg(eytr)) psitr=0.5d0*atan2(uutr,qqtr)*180.d0/pi chitr=0.5d0*asin(vvtr)*180.d0/pi ivac=2 anv=anvrfl xvrfl=xvout xv=xvout call vacuum_rt(xv,xvout,ivac) xv=xvout return 111 format(20(1x,e12.5)) end subroutine vacuum_rt(xvstart,xvend,ivac) use reflections, only : inters_linewall,inside implicit none integer*4 ivac integer nbb,nlim,i,imax parameter(nbb=5000) real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6) real*8 xv0(3) real*8 rlim(nbb),zlim(nbb) common/wrefl/walln common/limiter/rlim,zlim,nlim common/anv/anv common/dsds/dst common/psinv/psinv common/densbnd/psdbnd common/dstvac/dstvac c ivac=1 first interface plasma-vacuum c ivac=2 second interface vacuum-plasma after wall reflection c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection ! the real argument associated to xvstart is in a common block ! used by fwork!!! xv0=xvstart call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim), . nlim,smax,walln) smax=smax*1.d2 rrm=1.d-2*sqrt(xv0(1)**2+xv0(2)**2) zzm=1.d-2*xv0(3) if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then ! first wall interface is outside-inside if (dot_product(walln,walln)