SUBROUTINE AHWFNC WFNC001
c modified by leo mcnamara to handle no-field case as well
c
c DO NOT FORGET to set the magnetic pole at the geographic pole
c if you are using the no-field case ---- this is usually done
c automatically in raytrace_dp.for
c
c----- last modified 20-mar-87
implicit real*8 (a-h,o-z)
character*8 id,ndate,modx,mody,modz,modrin,koll,modsav
C CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE WFNC002
C APPLETON-HARTREE FORMULA WITH FIELD, NO COLLISIONS WFNC003
COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN WFNC004
COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,KAY2I, WFNC005
1 H,HI,PHPT,PHPTI,PHPR,PHPRI,PHPTH,PHPTHI,PHPPH,PHPPHIWFNC006
2 ,PHPOM,PHPOMI,PHPKR,PHPKRI,PHPKTH,PHPKTI,PHPKPH,PHPKPIWFNC007
3 ,KPHPK,KPHPKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN WFNC008
COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX WFNC009
COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRWFNC010
1 ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP WFNC011
COMMON /ZZ/ MODZ,Z(4) WFNC012
COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART WFNC013
COMMON R,TH,PH,KR,KTH,KPH /WW/ ID(10),W0,W(400) WFNC014
EQUIVALENCE (RAY,W(1)),(F,W(6)) WFNC015
LOGICAL SPACE WFNC016
REAL*8 KR,KTH,KPH,K2,KPHPK,KPHPKI,KAY2,KAY2I,N2,NNP,LPOLAR,LPOLRI WFNC017
real*8 k,logten
c
data ipass / 0 /
data u / 1.0d0 /
c
ENTRY RINDEX WFNC027
if(ray.eq.0.d0.and.ipass.eq.0) print*,' no magnetic field'
ipass = 1
OM=PIT2*1.d6*F WFNC028
C2=C*C WFNC029
K2=KR*KR+KTH*KTH+KPH*KPH WFNC030
OM2=OM*OM WFNC031
VR =C/OM*KR WFNC032
VTH=C/OM*KTH WFNC033
VPH=C/OM*KPH WFNC034
CALL ELECTX WFNC035
c
c allow for ray penetrating grid - 64 is chosen at random(!)
if(x.ge.64.d0) return
c
c treat no-field case separately
if(ray.ne.0.d0) go to 10
nnp = 1.
pnpx = -0.5d0
polari = 1.d0
lpolri = 1.d0
n2 = 1.d0 - x
cc type*,' ELECTX - h,x ', r - w(2),x
c if(x.lt.1.d-10) n2 = 1.d0
go to 20
10 CALL MAGY WFNC036
V2=VR**2+VTH**2+VPH**2 WFNC037
VDOTY=VR*YR+VTH*YTH+VPH*YPH WFNC038
YLV=VDOTY/V2 WFNC039
YL2=VDOTY**2/V2 WFNC040
YT2=Y**2-YL2 WFNC041
YT4=YT2*YT2 WFNC042
UX=U-X WFNC043
UX2=UX*UX WFNC044
RAD=RAY*dSQRT(YT4+4.d0*YL2*UX2) WFNC045
D=2.*UX-YT2+RAD WFNC046
D2=D*D WFNC047
N2=1.d0-2.d0*X*UX/D WFNC048
c
c avoid SQRT error in calculation of KAY2
cc if(n2.lt.0.d0) print*,' AHWFNC x,ux,d,n2 ',x,ux,d,n2
if(n2.lt.0.d0) n2 = 0.d0
c
PNPPS=2.d0*X*UX*(-1.d0+(YT2-2.d0*UX2)/RAD)/D2 WFNC049
PPSPR= YL2/Y*PYPR -(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV WFNC050
PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV WFNC051
PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV WFNC052
PNPX=-(2.d0*UX2-YT2*(U-2.d0*X)+(YT4*(U-2.d0*X)+4.d0*
a YL2*UX*UX2)/RAD)/D2 WFNC053
PNPY=2.d0*X*UX*(-YT2+(YT4+2.d0*YL2*UX2)/RAD)/(D2*Y) WFNC054
NNP=N2-(2.d0*X*PNPX+Y*PNPY) WFNC055
20 continue
PNPR =PNPX*PXPR
PNPTH=PNPX*PXPTH
PNPPH=PNPX*PXPPH
c bypass for no-field case
if(ray.eq.0.d0) go to 30
PNPR =PNPr +PNPY*PYPR +PNPPS*PPSPR WFNC056
PNPTH=PNPth +PNPY*PYPTH+PNPPS*PPSPTH WFNC057
PNPPH=PNPph +PNPY*PYPPH+PNPPS*PPSPPH WFNC058
PNPVR =PNPPS*(VR *YL2-VDOTY*YR )/V2 WFNC059
PNPVTH=PNPPS*(VTH*YL2-VDOTY*YTH)/V2 WFNC060
PNPVPH=PNPPS*(VPH*YL2-VDOTY*YPH)/V2 WFNC061
POLARI=dSQRT(V2)*(YT2-RAD)/(2.*VDOTY*UX) WFNC064
GAM=(-YT2+RAD)/(2.d0*UX) WFNC065
LPOLRI=X*dSQRT(YT2)/(UX*(U+GAM)) WFNC066
30 continue
PNPT=PNPX*PXPT WFNC062
SPACE=N2.EQ.1.d0 WFNC063
KAY2=OM2/C2*N2 WFNC067
IF(RSTART.EQ.0.d0) GO TO 1 WFNC068
c
c need to avoid SQRT error with 'scale'
scale = kay2/k2
if(scale.lt.0.d0) print*,' AHWFNC kay2,k2,scale',kay2,k2,scale
c if N2 = 0, scale = 0, & get 1/0 at WFNC039
if(scale.eq.0.d0) scale = 1.d0
if(scale.lt.0.d0) scale = 1.d0
c
SCALE=dSQRT(scale) WFNC069
KR =SCALE*KR WFNC070
KTH=SCALE*KTH WFNC071
KPH=SCALE*KPH WFNC072
1 CONTINUE WFNC073
C********* CALCULATES A HAMILTONIAN H WFNC074
c H=.5d0*(C2*K2/OM2-N2) WFNC075
h = k2 / om2
h = h * c2
h = 0.5d0 * ( h - n2 )
C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO WFNC076
C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA, AND KPHI. WFNC077
PHPT =-PNPT WFNC078
PHPR =-PNPR WFNC079
c???? phpr = - pnpr * dsqrt(n2)
c???? check this some time Leo 26-may-87
PHPTH=-PNPTH WFNC080
c???? phpth = - pnpth * dsqrt(n2)
PHPPH=-PNPPH WFNC081
PHPOM=-NNP/OM WFNC082
PHPKR =C2/OM2*KR
PHPKTH=C2/OM2*KTH
PHPKPH=C2/OM2*KPH
if(ray.eq.0.d0) go to 40
PHPKR =phpKR -C/OM*PNPVR WFNC083
PHPKTH=phpKTH -C/OM*PNPVTH WFNC084
PHPKPH=phpKPH -C/OM*PNPVPH WFNC085
40 continue
KPHPK=N2 WFNC086
RETURN WFNC087
END WFNC088-
block data appleton
character*8 id,ndate,modx,mody,modz,modrin,koll,modsav
COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN WFNC004
COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,KAY2I, WFNC005
1 H,HI,PHPT,PHPTI,PHPR,PHPRI,PHPTH,PHPTHI,PHPPH,PHPPHIWFNC006
2 ,PHPOM,PHPOMI,PHPKR,PHPKRI,PHPKTH,PHPKTI,PHPKPH,PHPKPIWFNC007
3 ,KPHPK,KPHPKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN WFNC008
COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX WFNC009
COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRWFNC010
1 ,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP WFNC011
COMMON /ZZ/ MODZ,Z(4) WFNC012
COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART WFNC013
COMMON R,TH,PH,KR,KTH,KPH /WW/ ID(10),W0,W(400) WFNC014
EQUIVALENCE (RAY,W(1)),(F,W(6)) WFNC015
LOGICAL SPACE WFNC016
REAL*8 KR,KTH,KPH,K2,KPHPK,KPHPKI,KAY2,KAY2I,N2,NNP,LPOLAR,LPOLRI WFNC017
real*8 k,logten
data coll,kay2i,hi,phpti,phpri,phpthi,phpphi,phpomi,phpkri,
a phpkti,phpkpi,kphpki,polar,lpolar,x,pxpr,pxpth,pxpph,pxpt,
a y,pypr,pypth,pypph,yr,pyrpr,pyrpt,pyrpp,yth,pytpr,pytpt,
a pytpp,yph,pyppr,pyppt,pyppp / 35*0.d0 /
data field / 1.d0 /
data modrin,modz / 'appleton','-hartree',' formula',' ' /
end