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