From 5f8331937dbe53821ee58f017404baac2e7ad0b9 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 18 Apr 2013 15:01:41 +0000 Subject: [PATCH 1/8] Branch for OXB analysis created From ea75b096cb352a6abe0a0fd792720dda93db32ee Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 18 Apr 2013 15:26:27 +0000 Subject: [PATCH 2/8] Added prints for OXB --- src/gray.f | 81 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/src/gray.f b/src/gray.f index 46505e2..c8ab213 100644 --- a/src/gray.f +++ b/src/gray.f @@ -153,6 +153,7 @@ c if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES' if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm + call print_tox end if c c compute power and current density profiles for all rays @@ -267,6 +268,7 @@ 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 + call calc_tox(i) if(ipolc.eq.0.and.iov(j,k).eq.1) then call pol_limit(qqin,uuin,vvin) ipolc=1 @@ -355,6 +357,77 @@ c end c c +c + subroutine calc_tox(i) + implicit real*8 (a-h,o-z) +c computation of O-X transmission function +c calculated only for xg=(wpl/w)**2>0.5 +c values stored at maximum xg and at maximum transmission + parameter(pi=3.14159265358979d0) + dimension derxg(3),deryg(3),xv(3),anv(3) +c + common/xgxg/xg + common/ygyg/yg + common/xv/xv + common/anv/anv + common/nplr/anpl,anpr + common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/nprw/anprr,anpri + common/parwv/ak0,akinv,fhz + common/dxgyg/derxg,deryg + common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx + common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx + + anplopt=sqrt(yg/(yg+1.0d0)) + thkb=1.8d2/pi*acos(anpl/sqrt(an2)) + dlen=sqrt(derxg(1)**2+derxg(2)**2+derxg(3)**2) + if (dlen.gt.0.0d0) then + dlen=xg/dlen + tox=exp(-pi*ak0*dlen*sqrt(0.5d0*yg)* + . (2*(1.0d0+yg)*(anplopt-abs(anpl))**2+anpr**2)) + else + tox=0.0d0 + end if + if (xg.gt.0.5d0.and.xg.le.1.0d0) then + if (xg.gt.xgmax) then + ixmx=i + toxxmx=tox + xgmax=xg + ygxmx=yg + anplxmx=anpl + anprxmx=anpr + thkbxmx=thkb + end if + if (tox.gt.toxmax) then + itmx=i + toxmax=tox + xgtmx=xg + ygtmx=yg + anpltmx=anpl + anprtmx=anpr + thkbtmx=thkb + end if + write(46,98) i,tox,xg,yg,anpl,anpr,thkb,xv/1.d2,anv,derxg,deryg + end if +c + return + 98 format(1x,i8,30(1x,e16.8e3)) + end +c +c +c + subroutine print_tox + implicit real*8 (a-h,o-z) + common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx + common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx + write(44,98) ixmx,toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx + write(45,98) itmx,toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx +c + return + 98 format(1x,i8,30(1x,e16.8e3)) + end +c +c c subroutine print_output(i,j,k) implicit real*8 (a-h,o-z) @@ -2585,6 +2658,9 @@ c common/ierr/ierr common/istop/istop common/ipol/ipolc +c used only for computation of O-X transmission function + common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx + common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx c istpr=0 istpl=1 @@ -2592,6 +2668,11 @@ c istep=0 istop=0 ipolc=0 +c + ixmx=0 + itmx=0 + xgmax=0.0d0 + toxmax=0.0d0 c return end From fbe86719f31374dd3ff9a96d22235934e4e5ff09 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 9 May 2013 09:18:53 +0000 Subject: [PATCH 3/8] local modifications in gaussfit copied in oxb version --- src/gray.f | 453 ++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 416 insertions(+), 37 deletions(-) diff --git a/src/gray.f b/src/gray.f index c8ab213..7f1ee7e 100644 --- a/src/gray.f +++ b/src/gray.f @@ -578,9 +578,11 @@ c .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' write(8,*) ' #istep j k xt yt zt rt psin' write(9,*) ' #istep j k xt yt zt rt psin' + write(82,*) ' #istep j k xt yt zwspl zwparab' write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' - write(12,*) ' #i sst psi w1 w2' + write(12,*) ' #i sst w1 w2 phiw rci1 rci2 phir errw errr '// + . 'dk1 dk2 dkpar phik dnpar' write(7,*)'#Icd Pa Jphimx dPdVmx '// .'rhotj rhotjava rhotp rhotpav '// .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' @@ -624,6 +626,7 @@ c common/warm/iwarm,ilarm common/ieccd/ieccd common/idst/idst + common/iplane/iplane c common/filesn/filenmeqq,filenmprf,filenmbm c @@ -665,6 +668,8 @@ c common/mode/sox common/angles/alpha0,beta0 common/scal/iscal + + common/waist/w0csi,w0eta c open(2,file='gray.data',status= 'unknown') c @@ -742,8 +747,9 @@ 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 iplane=0/1 1 plane integration wavefronts c - read(2,*) dst,nstep,istpr0,istpl0,idst + read(2,*) dst,nstep,istpr0,istpl0,idst,iplane read(2,*) igrad,ipass,rwallm read(2,*) iox, psipol0,chipol0 c @@ -2953,6 +2959,7 @@ c hh=h*0.5d0 h6=h/6.0d0 c + ttest=0 do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 @@ -2970,17 +2977,17 @@ c ddgr(iv,jv)=ggri(iv,jv,j,k) end do end do - call fwork(yy,fk2) + call fwork(j,k,1,yy,fk2) c do ieq=1,ndim yy(ieq)=y(ieq)+fk2(ieq)*hh end do - call fwork(yy,fk3) + call fwork(j,k,2,yy,fk3) c do ieq=1,ndim yy(ieq)=y(ieq)+fk3(ieq)*h end do - call fwork(yy,fk4) + call fwork(j,k,3,yy,fk4) c do ieq=1,ndim ywrk(ieq,j,k)=y(ieq) @@ -3034,7 +3041,7 @@ c end do end if c - call fwork(yy,yyp) + call fwork(j,k,3,yy,yyp) c do ieq=1,ndim ypwrk(ieq,j,k)=yyp(ieq) @@ -3047,14 +3054,18 @@ c c c c - subroutine fwork(y,dery) + subroutine fwork(j,k,isubst,y,dery) implicit real*8 (a-h,o-z) parameter(ndim=6) dimension y(ndim),dery(ndim) 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) + dimension yy11(6,3),derdxv11(3,3) + dimension dery0(3) c + save yy11,derdxv11 + common/gr/gr2 common/dgr/dgr2,dgr,ddgr common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 @@ -3072,6 +3083,10 @@ c common/anv/anv common/xv/xv common/idst/idst +c + common/iplane/iplane +c + common/dery0/dery0,dery0mod c xx=y(1) yy=y(2) @@ -3179,6 +3194,14 @@ c c derdnm=derdnm+derdnv(iv)**2 end do + + if(j.eq.1.and.k.eq.1) then + do ii=1,3 + yy11(ii,isubst)=y(ii) + yy11(ii+3,isubst)=y(ii+3) + derdxv11(ii,isubst)=derdxv(ii) + enddo + endif c derdnm=sqrt(derdnm) c @@ -3195,7 +3218,16 @@ 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)) + if (iplane.eq.1.and.j.gt.1) then +c advance outer rays to the plane through x_11 and perp to N_11 + denom=-((yy11(4,isubst)*derdnv(1)+yy11(5,isubst)*derdnv(2) + . +yy11(6,isubst)*derdnv(3)) + . -(derdxv11(1,isubst)*(y(1)-yy11(1,isubst))+ + . derdxv11(2,isubst)*(y(2)-yy11(2,isubst))+ + . derdxv11(3,isubst)*(y(3)-yy11(3,isubst)))) + else + denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3)) + end if end if c c coefficient for integration in s @@ -3208,6 +3240,14 @@ c dery(4) = derdxv(1)/denom dery(5) = derdxv(2)/denom dery(6) = derdxv(3)/denom + + if(j.eq.1) then + do ll=1,3 + dery0(ll)=dery(ll) + enddo + dery0mod=sqrt(dery0(1)*dery0(1)+ + . dery0(2)*dery0(2)+dery0(3)*dery0(3)) + endif c c vgv : ~ group velocity c @@ -5742,34 +5782,177 @@ c gg=F(u)/u with F(u) as in Cohen paper subroutine projxyzt(iproj,nfile) - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36) + implicit real*8 (a-h,o-z) + parameter(jmx=31,kmx=36,nmx=8000) + parameter(pi=3.14159265358979d0) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) + complex*16 ui,aac,bbc,ccc,ddc,aak,bbk,cck + parameter(ui=(0.0d0,1.0d0)) + dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) + parameter(nrmax=(jmx-1)*kmx+1) + dimension xtiv(nrmax),ytiv(nrmax),zwjv(nrmax),w(nrmax) + + dimension pvett(3),pvettn(3),dery0(3),dery0n(3) + dimension avn(3,jmx,kmx) + dimension aplane(3,jmx,kmx),asip(jmx,kmx) + dimension gri(3,jmx,kmx) + +c parameter(nxmax=2*jmx-1) + parameter(nxmax=2*kmx) + parameter(kspl=3,nxest=nxmax+4) + parameter(km=kspl+1,nu=nxest-km) + parameter(ne=nxest,kb1=kspl*nu+km,kb2=kb1+nu-kspl) + parameter(lwrk1=nu*nu*(2+kb1+kb2)+ + . 2*(2*nu+km*(nrmax+ne)+ne-2*kspl)+kb2+1) + parameter(lwrk2=nu*nu*(kb2+1)+kb2) + parameter(kwrk=nrmax+(nxest-kspl-km)*(nxest-kspl-km)) + dimension xgridv(nxmax),ygridv(nxmax) + dimension zwint(nxmax*nxmax),ccexp(nxmax*nxmax) + dimension tx(nxest),ty(nxest),wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk) + parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax) + dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) c common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk - common/psinv11/psinv11 + common/derip11/dpsi11dr,dpsi11dz,psinv11,bvx11,bvy11,bvz11 common/istep/istep + common/rwmax/rwmax common/ss/st + common/dnpar/dnpara + common/atjki/tauv,alphav + common/waist/w0csi,w0eta + + common/dery0/dery0,dery0mod + common/iplane/iplane + common/gradjk/gri c - rtimn=1.d+30 - rtimx=-1.d-30 + x4m=0.0d0 + x3ym=0.0d0 + x2y2m=0.0d0 + xy3m=0.0d0 + y4m=0.0d0 + x2zrm=0.0d0 + xyzrm=0.0d0 + y2zrm=0.0d0 + x2zwm=0.0d0 + xyzwm=0.0d0 + y2zwm=0.0d0 + z2wm=0.0d0 + z2rm=0.0d0 +c initialize grid dimension for spline interpolation + xmaxgrid=2*max(w0csi,w0eta) + iray=0 + if(iplane.le.1) then + do j=1,nrayr + kktx=nrayth + if(j.eq.1) kktx=1 + do k=1,kktx + do ll=1,3 + aplane(ll,j,k)=ywrk(ll,j,k) + enddo + asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2 + enddo + enddo + endif + if(iplane.eq.2) then +c prjection parallel to vg on the plane perpendicular to n0 passing through x11 + do j=2,nrayr + do k=1,nrayth + avmod=sqrt(ypwrk(1,j,k)**2+ypwrk(2,j,k)**2+ypwrk(3,j,k)**2) + do ii=1,3 + avn(ii,j,k)=ypwrk(ii,j,k)/avmod + enddo + deltaxdotn0=(ywrk(1,1,1)-ywrk(1,j,k))*ywrk(4,1,1)+ + . (ywrk(2,1,1)-ywrk(2,j,k))*ywrk(5,1,1)+ + . (ywrk(3,1,1)-ywrk(3,j,k))*ywrk(6,1,1) + avndotn0=avn(1,j,k)*ywrk(4,1,1)+ + . avn(2,j,k)*ywrk(5,1,1)+ + . avn(3,j,k)*ywrk(6,1,1) + aalpha=deltaxdotn0/avndotn0 + do ll=1,3 + aplane(ll,j,k)=ywrk(ll,j,k)+aalpha*avn(ll,j,k) + enddo + enddo + enddo + do ll=1,3 + aplane(ll,1,1)=ywrk(ll,1,1) + enddo + do j=1,nrayr + kktx=nrayth + if(j.eq.1) kktx=1 + do k=1,kktx + asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2 + enddo + enddo + endif + if(iplane.eq.3) then +c ortogonal projection on the plane perpendicular to n0 passing through x11 + do j=2,nrayr + do k=1,nrayth + deltaxdotn0=(ywrk(1,1,1)-ywrk(1,j,k))*ywrk(4,1,1)+ + . (ywrk(2,1,1)-ywrk(2,j,k))*ywrk(5,1,1)+ + . (ywrk(3,1,1)-ywrk(3,j,k))*ywrk(6,1,1) + an02=ywrk(4,1,1)**2+ywrk(5,1,1)**2+ywrk(6,1,1)**2 + aalpha=deltaxdotn0/an02 + do ll=1,3 + aplane(ll,j,k)=ywrk(ll,j,k)+aalpha*ywrk(ll+3,1,1) + enddo + enddo + enddo + do ll=1,3 + aplane(ll,1,1)=ywrk(ll,1,1) + enddo +c Si evaluation on the projection plane(Taylor, first order) + do j=1,nrayr + kktx=nrayth + if(j.eq.1) kktx=1 + do k=1,kktx + asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2+ + . gri(1,j,k)*(aplane(1,j,k)-ywrk(1,j,k))+ + . gri(2,j,k)*(aplane(2,j,k)-ywrk(2,j,k))+ + . gri(3,j,k)*(aplane(3,j,k)-ywrk(3,j,k)) + enddo + enddo + endif +c + do j=1,nrayr + kktx=nrayth + if(j.eq.1) kktx=1 + do k=1,kktx + zwj=asip(j,k)+ + . 0.5*(tauv(j,k,istep)-tauv(1,1,istep)) 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) + dx=aplane(1,j,k)-aplane(1,1,1) + dy=aplane(2,j,k)-aplane(2,1,1) + dz=aplane(3,j,k)-aplane(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) + + if (j>1) then + k2=mod(k+kktx/4-1,kktx)+1 + dx2=aplane(1,j,k2)-aplane(1,1,1) + dy2=aplane(2,j,k2)-aplane(2,1,1) + dz2=aplane(3,j,k2)-aplane(3,1,1) + pvett(1)=dy*dz2-dy2*dz + pvett(2)=dz*dx2-dz2*dx + pvett(3)=dx*dy2-dx2*dy + pvettmod=sqrt(pvett(1)**2+pvett(2)**2+pvett(3)**2) + do ll=1,3 + pvettn(ll)=pvett(ll)/pvettmod + dery0n(ll)=dery0(ll)/dery0mod + enddo +c write(*,*) 'dotn0',j,k,(pvettn(1)*dirx+ +c . pvettn(2)*diry+pvettn(3)*dirz)/dir +c write(*,*) 'dotvg0',j,k,(pvettn(1)*dery0(1)+ +c . pvettn(2)*dery0(2)+pvettn(3)*dery0(3))/dery0mod + endif +c dirx=dery0(1) +c diry=dery0(2) +c dirz=dery0(3) +c dir=dery0mod c if(j.eq.1.and.k.eq.1) then csth1=dirz/dir @@ -5784,39 +5967,235 @@ c 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) + rti=sqrt(xti**2+yti**2) +c store x,y,z values for spline interpolation + iray=iray+1 + xtiv(iray)=xti + ytiv(iray)=yti + zwjv(iray)=zwj +c initialize grid dimension for spline interpolation + xmaxgrid=max(xmaxgrid,rti) 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 + bxti= bvx11*csps1-bvy11*snps1 + byti=(bvx11*snps1+bvy11*csps1)*csth1-bvz11*snth1 +c bzti=(bvx11*snps1+bvy11*csps1)*snth1+bvz11*csth1 +c + rr11=sqrt(ywrk(1,1,1)**2+ywrk(2,1,1)**2) + dpsx=dpsi11dr*ywrk(1,1,1)/rr11 + dpsy=dpsi11dr*ywrk(2,1,1)/rr11 + dpsz=dpsi11dz + xdpsi=dpsx*csps1-dpsy*snps1 + ydpsi=(dpsx*snps1+dpsy*csps1)*csth1-dpsz*snth1 + zdpsi=(dpsx*snps1+dpsy*csps1)*snth1+dpsz*csth1 + sq=sqrt(xdpsi**2+ydpsi**2+zdpsi**2) + if(sq.gt.0.0d0) then + xdpsi=dx*xdpsi/sq + ydpsi=dx*ydpsi/sq + zdpsi=dx*zdpsi/sq + end if c if(k.eq.1) then xti1=xti yti1=yti zti1=zti rti1=rti - end if + xdpsi1=xdpsi + ydpsi1=ydpsi + zdpsi1=zdpsi + end if +c + if(istep.eq.0) + . write(10,111) istep,j,k,xti,yti,zti,dirxt,diryt,dirzt,dir +c dr=sqrt(dx**2+dy**2+dz**2) +c write(11,111) istep,j,k,dx/dr,dy/dr,dz/dr, +c . snth1*snps1,snth1*csps1,csth1, +c . snth1*snps1*dx/dr+snth1*csps1*dy/dr+csth1*dz/dr - 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 + if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nrayr)) + . write(nfile,111) istep,j,k,xti,yti,zti,rti, + . xdpsi,ydpsi,zdpsi c + x4m=x4m+xti**4 + x3ym=x3ym+xti**3*yti + x2y2m=x2y2m+xti**2*yti**2 + xy3m=xy3m+xti*yti**3 + y4m=y4m+yti**4 + x2zrm=x2zrm+xti**2*zti + xyzrm=xyzrm+xti*yti*zti + y2zrm=y2zrm+yti**2*zti + x2zwm=x2zwm+xti**2*zwj + xyzwm=xyzwm+xti*yti*zwj + y2zwm=y2zwm+yti**2*zwj + z2wm=z2wm+zwj*zwj + z2rm=z2rm+zti*zti 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 +c + if((iproj.eq.1.and.j.gt.1).or.(iproj.eq.0.and.j.eq.nrayr)) + . write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1, + . xdpsi1,ydpsi1,zdpsi1 if(iproj.eq.1) write(nfile,*) ' ' end do -c +c write(nfile,*) ' ' c - write(12,99) istep,st,psinv11,rtimn,rtimx +c computation of the SI paraboloid +c + denw= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - + . x2y2m*(2*x3ym*xy3m + x4m*y4m) + aaw= -(x2y2m*xy3m*xyzwm) + x2y2m**2*y2zwm - + . x3ym*xy3m*y2zwm + x3ym*xyzwm*y4m + + . x2zwm*(xy3m**2 - x2y2m*y4m) + ccw= x2y2m**2*xyzwm + x4m*xy3m*y2zwm - + . x2y2m*(x2zwm*xy3m + x3ym*y2zwm) + x2zwm*x3ym*y4m - + . x4m*xyzwm*y4m + bbw= x2y2m**2*x2zwm - x2zwm*x3ym*xy3m + x4m*xy3m*xyzwm + + . x3ym**2*y2zwm - x2y2m*(x3ym*xyzwm + x4m*y2zwm) + aaw=aaw/denw + bbw=bbw/denw + ccw=ccw/denw + phiw = 0.5d0*atan(ccw/(aaw-bbw+1.0d-32)) + phiwdeg=phiw*180.d0/pi + aaaw = 0.5d0*(aaw+bbw + (aaw-bbw)/cos(2.0d0*phiw)) + bbbw = aaw+bbw - aaaw + errw= 2.0d0*aaw*bbw*x2y2m + ccw**2*x2y2m - 2.0d0*aaw*x2zwm + + - 2.0d0*aaw*ccw*x3ym + aaw**2*x4m + 2.0d0*bbw*ccw*xy3m - + - 2.0d0*ccw*xyzwm - 2.0d0*bbw*y2zwm + bbw**2*y4m + z2wm + errw=sqrt(abs(errw)/dble((nrayr-1)*nrayth)) + wcsi = 1.0d0/sqrt(aaaw) + weta = 1.0d0/sqrt(bbbw) +c +c computation of paraboloid Sr=z-(ax^2+cxy+by^2)=const +c + denr= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - + . x2y2m*(2*x3ym*xy3m + x4m*y4m) + aar= -(x2y2m*xy3m*xyzrm) + x2y2m**2*y2zrm - + . x3ym*xy3m*y2zrm + x3ym*xyzrm*y4m + + . x2zrm*(xy3m**2 - x2y2m*y4m) + ccr= x2y2m**2*xyzrm + x4m*xy3m*y2zrm - + . x2y2m*(x2zrm*xy3m + x3ym*y2zrm) + x2zrm*x3ym*y4m - + . x4m*xyzrm*y4m + bbr= x2y2m**2*x2zrm - x2zrm*x3ym*xy3m + x4m*xy3m*xyzrm + + . x3ym**2*y2zrm - x2y2m*(x3ym*xyzrm + x4m*y2zrm) + aar=aar/denr + bbr=bbr/denr + ccr=ccr/denr + phir = 0.5d0*atan(ccr/(aar-bbr+1.0d-32)) + phirdeg=phir*180.d0/pi + aaar = 0.5d0*(aar+bbr + (aar-bbr)/cos(2.0d0*phir)) + bbbr = aar+bbr - aaar + errr= 2.0d0*aar*bbr*x2y2m + ccr**2*x2y2m - 2.0d0*aar*x2zrm + + - 2.0d0*aar*ccr*x3ym + aar**2*x4m + 2.0d0*bbr*ccr*xy3m - + - 2.0d0*ccr*xyzrm - 2.0d0*bbr*y2zrm + bbr**2*y4m + z2rm + errr=sqrt(abs(errr)/dble((nrayr-1)*nrayth)) + rcicsi=-2.0d0*aaar + rcieta=-2.0d0*bbbr +c +c computation of Fourier Transform of exp[i k0 (Sr + i Si)] +c k spectrum +c + aar=-ak0*aar + bbr=-ak0*bbr + ccr=-ak0*ccr + aac=aar+ui*aaw + bbc=bbr+ui*bbw + ccc=ccr+ui*ccw + ddc=ccc*ccc-4.0d0*aac*bbc + aak=bbc/ddc + bbk=aac/ddc + cck=-ccc/ddc + akw=dimag(aak) + bkw=dimag(bbk) + ckw=dimag(cck) + dkpar2=4.0d0*(bxti**2*bkw+byti**2*akw-bxti*byti*ckw)/ + . (4.0d0*akw*bkw-ckw*ckw) + phik = 0.5d0*atan(ckw/(akw-bkw+1.0d-32)) + aakw = 0.5d0*(akw+bkw + (akw-bkw)/cos(2.0d0*phik)) + bbkw = akw+bkw - aakw + dk1 = 1.0d0/sqrt(aakw) + dk2 = 1.0d0/sqrt(bbkw) +c co = cos(phik) +c si = sin(phik) +c bxtir= co*bxti+si*byti +c bytir=-si*bxti+co*byti +c dkpar2=bxtir**2/aakw+bytir**2/bbkw + dkpar=sqrt(dkpar2) + dnpar=dkpar/ak0 + dnpara=0.5d0*dnpar +c +c spectral distribution of electric field E(k): +c exp[-(N//-N//0)^2/(DNPAR^2)] +c dnpar^2 = 2 (Delta N//)^2 , Maj !!! +c +c in vacuum dkpar^2 = 4/w0^2 +c +c in common dnpara=dnpar/2 to match westerhof Delta function +c Delta -> exp[-(...)^2/(2*DeltaQ)] +c + + +c Spline fit + if(nrayr.gt.1) then + + npoints=iray + xmin=-xmaxgrid + xmax=xmaxgrid + ymin=-xmaxgrid + ymax=xmaxgrid + nxgrid=2*nrayr-1 + dxgrid=(xmax-xmin)/(nxgrid-1) + do i=1,nxgrid + xgridv(i)=xmin+(i-1)*dxgrid + ygridv(i)=xgridv(i) + end do + +c call interp spline + iopt=0 + sspl=1.0d-3 + eps=1.0d-7 + do iray=1,npoints + w(iray) = 1.0d0 + end do + call surfit(iopt,npoints,xtiv,ytiv,zwjv,w,xmin,xmax,ymin,ymax, + . kspl,kspl,sspl,nxest,nxest,nxest,eps, + . nkntx,tx,nknty,ty,ccexp,resid,wrk1,lwrk1,wrk2,lwrk2, + . iwrk,kwrk,ierr) + if (ierr.gt.0) then + print*, 'surfit:',istep,nkntx,nknty,ierr,resid + do j=1,nxgrid + do i=1,nxgrid + zwint(nxgrid*(i-1)+j)=0.0d0 + end do + end do + else +c call eval spline + call bispev(tx,nkntx,ty,nknty,ccexp,kspl,kspl,xgridv,nxgrid, + . ygridv,nxgrid,zwint,wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ierr) + if (ierr.ne.0) print*, 'bispev:',istep,ierr + end if + do j=1,nxgrid + do i=1,nxgrid + ij=nxgrid*(i-1)+j + write(82,111) istep,i,j,xgridv(i),ygridv(j),zwint(ij), + . aaw*xgridv(i)**2+ccw*xgridv(i)*ygridv(j)+bbw*ygridv(j)**2 + end do + write(82,*) '' + end do + write(82,*) '' + end if + + write(12,99) istep,st, + . wcsi,weta,phiwdeg,rcicsi,rcieta,phirdeg,errw,errr, + . dk1,dk2,dkpar,phik*180.0d0/pi,dnpar +c return - 99 format(i5,12(1x,e16.8e3)) + 99 format(i5,22(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) end + c c c @@ -6561,7 +6940,7 @@ c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection y(4)=anv(1) y(5)=anv(2) y(6)=anv(3) - call fwork(y,dery) + call fwork(1,1,3,y,dery) if (st.ge.smax.or.(psinv.gt.0.0d0.and.psinv.lt.psdbnd)) exit i=i+1 end do From e74ec6eecbe07f9667e3a0bcc863de6e4c8fd1b2 Mon Sep 17 00:00:00 2001 From: Alberto Mariani Date: Tue, 21 May 2013 15:41:15 +0000 Subject: [PATCH 4/8] projxyzt modified:in output correct inverse radii of curvature --- src/gray.f | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/src/gray.f b/src/gray.f index 16e3505..f23c239 100644 --- a/src/gray.f +++ b/src/gray.f @@ -503,8 +503,8 @@ c write(4,*)' #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(8,*) ' #istep j k xt yt zt rt psin' - write(9,*) ' #istep j k xt yt zt rt psin' + write(8,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn' + write(9,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn' write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' write(12,*) ' #i sst psi w1 w2' @@ -2855,6 +2855,7 @@ c dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + dimension vgv(3),vgv11(3) c common/nray/nrayr,nrayth common/dsds/dst @@ -2867,6 +2868,9 @@ c common/gr/gr2 common/dgr/dgr2,dgr,ddgr common/igrad/igrad + + + common/vgv11/vgv,vgv11 c h=dst hh=h*0.5d0 @@ -2906,6 +2910,7 @@ c . +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq)) end do end do + if(j.eq.1) vgv11=vgv end do c call updatepos @@ -2970,9 +2975,10 @@ c implicit real*8 (a-h,o-z) parameter(ndim=6) dimension y(ndim),dery(ndim) - dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3) + dimension xv(3),anv(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) + dimension vgv(3),vgv11(3) c common/gr/gr2 common/dgr/dgr2,dgr,ddgr @@ -2991,6 +2997,7 @@ c common/anv/anv common/xv/xv common/idst/idst + common/vgv11/vgv,vgv11 c xx=y(1) yy=y(2) @@ -5655,12 +5662,14 @@ c gg=F(u)/u with F(u) as in Cohen paper implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) + dimension vgv(3),vgv11(3) c common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk common/psinv11/psinv11 common/istep/istep common/ss/st + common/vgv11/vgv,vgv11 c rtimn=1.d+30 rtimx=-1.d-30 @@ -5676,10 +5685,16 @@ c 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) + dirxn=ywrk(4,j,k) + diryn=ywrk(5,j,k) + dirzn=ywrk(6,j,k) + dirn=sqrt(dirxn*dirxn+diryn*diryn+dirzn*dirzn) + dirx=vgv11(1) + diry=vgv11(2) + dirz=vgv11(3) dir=sqrt(dirx*dirx+diry*diry+dirz*dirz) +c + vgdotn=dirxn*dirx+diryn*diry+dirzn*dirz c if(j.eq.1.and.k.eq.1) then csth1=dirz/dir @@ -5708,7 +5723,7 @@ c end if if(.not.(iproj.eq.0.and.j.eq.1)) - . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 + . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11,dir,dirn,vgdotn c if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti From 81dc2d0a3c03da8ae77b4ebe75bca406b431c228 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Fri, 20 Sep 2013 14:48:53 +0000 Subject: [PATCH 5/8] added print og psi,X,Y on R,z grid --- src/gray.f | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/src/gray.f b/src/gray.f index 7f1ee7e..47dd34f 100644 --- a/src/gray.f +++ b/src/gray.f @@ -910,7 +910,11 @@ c versus psi, rhop, rhot end if if (iequil.eq.1) call surf_anal - +c +c print psi,X,Y on the R,z grid +c + call print_2d +c if (iequil.ne.2.or.ipass.lt.0) then c set simple limiter as two cylindrical walls at rwallm and r00 nlim=5 @@ -1642,6 +1646,34 @@ c end c c +c + subroutine print_2d + implicit real*8 (a-h,o-z) + parameter(nnw=501,nnh=501) + dimension rv(nnw),zv(nnh),psin(nnw,nnh) + dimension btotal(nnw,nnh) +c + common/cpsin/rv,zv,psin + common/psival/psinv + common/eqnn/nr,nz,npp,nintp + common/parpl/brr,bphi,bzz,ajphi + common/xgxg/xg + common/ygyg/yg +c + do j=1,nr + rjcm=rv(j)*1.d2 + do k=1,nz + zkcm=zv(k)*1.d2 + call plas_deriv(rjcm,0.d0,zkcm) + write(91,111) rv(j),zv(k),psinv,brr,bphi,bzz,xg,yg + enddo + write(91,*) ' ' + enddo + return + 111 format(12(1x,e12.5)) + end +c +c c subroutine print_prof implicit real*8 (a-h,o-z) From 85b13a4b8d9dd1ed00cbcfce3031281bb39cd22f Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Fri, 11 Oct 2013 15:10:56 +0000 Subject: [PATCH 6/8] cleaned output. merged correction from trunk rev 52 (spline coefficients dimensions) --- src/gray.f | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/src/gray.f b/src/gray.f index 47dd34f..9883f15 100644 --- a/src/gray.f +++ b/src/gray.f @@ -376,7 +376,7 @@ c common/parwv/ak0,akinv,fhz common/dxgyg/derxg,deryg common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx - common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx +! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx anplopt=sqrt(yg/(yg+1.0d0)) thkb=1.8d2/pi*acos(anpl/sqrt(an2)) @@ -398,16 +398,16 @@ c anprxmx=anpr thkbxmx=thkb end if - if (tox.gt.toxmax) then - itmx=i - toxmax=tox - xgtmx=xg - ygtmx=yg - anpltmx=anpl - anprtmx=anpr - thkbtmx=thkb - end if - write(46,98) i,tox,xg,yg,anpl,anpr,thkb,xv/1.d2,anv,derxg,deryg +! if (tox.gt.toxmax) then +! itmx=i +! toxmax=tox +! xgtmx=xg +! ygtmx=yg +! anpltmx=anpl +! anprtmx=anpr +! thkbtmx=thkb +! end if +! write(46,98) i,tox,xg,yg,anpl,anpr,thkb,xv/1.d2,anv,derxg,deryg end if c return @@ -419,9 +419,9 @@ c subroutine print_tox implicit real*8 (a-h,o-z) common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx - common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx +! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx write(44,98) ixmx,toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx - write(45,98) itmx,toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx +! write(45,98) itmx,toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx c return 98 format(1x,i8,30(1x,e16.8e3)) @@ -586,6 +586,7 @@ c write(7,*)'#Icd Pa Jphimx dPdVmx '// .'rhotj rhotjava rhotp rhotpav '// .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' + write(44,*) '#i Tox X Y Npl Nperp thetakB' write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' else @@ -2698,7 +2699,7 @@ c common/ipol/ipolc c used only for computation of O-X transmission function common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx - common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx +! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx c istpr=0 istpl=1 @@ -2708,9 +2709,9 @@ c ipolc=0 c ixmx=0 - itmx=0 +! itmx=0 xgmax=0.0d0 - toxmax=0.0d0 +! toxmax=0.0d0 c return end @@ -3555,7 +3556,7 @@ c 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(lw10),cc10(lw01),cc02(lw02),cc20(lw20),cc11(lw11) + dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) dimension tfp(nrest),cfp(nrest),wrkfd(nrest) c common/eqnn/nr,nz,npp,nintp From 5a5406a5cf831eeb4480b4344f44b6a4c0430a15 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Tue, 10 Jun 2014 14:03:41 +0000 Subject: [PATCH 7/8] fix to initialisation in reflections module propagated to other branches --- src/reflections.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/reflections.f90 b/src/reflections.f90 index ffe59bd..6e04214 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -32,7 +32,8 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) integer :: i,j,ni,iint real(r8), dimension(2) :: si,ti real(r8) :: drw,dzw,xint,yint,rint,l,kxy - real(r8) :: tol=sqrt(epsilon(1.0_r8)) + real(r8) :: tol + tol=sqrt(epsilon(1.0_r8)) sint=huge(sint) iint=0 normw=0.0_r8 From 4909c9fe13ab78b4f3578e4f54b14f41953a3887 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Tue, 10 Jun 2014 14:03:41 +0000 Subject: [PATCH 8/8] fix to initialisation in reflections module propagated to other branches --- src/reflections.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/reflections.f90 b/src/reflections.f90 index ffe59bd..6e04214 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -32,7 +32,8 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) integer :: i,j,ni,iint real(r8), dimension(2) :: si,ti real(r8) :: drw,dzw,xint,yint,rint,l,kxy - real(r8) :: tol=sqrt(epsilon(1.0_r8)) + real(r8) :: tol + tol=sqrt(epsilon(1.0_r8)) sint=huge(sint) iint=0 normw=0.0_r8