C THIS ROUTINE WILL BE USED FOR JUPITER RAY TRACING MODIFIED 9/20/94 C THIS VERSION DIFFERS FROM "JUPITER" BECAUSE IT CONTAINS PARAMETERS C TO MODIFY THE IO TORUS. C C THIS ROUTINE CONTAINS THE DIVINE AND GARRETT PLASMA MODEL C AND THE NEW O6 MAGNETIC FIELD MODEL C C IT ALSO CONTAINS THE NEW FCT ROUTINE WHICH IS THE SAME AS IN C "NRAY.FOR", I.E. IT REDEFINES THE EXPRESSIONS IN ORDER TO ELIMINATE C THE ARTIFICIAL SINGULARITY AT F=FG. C C C THE INPUT FOR STARTING POSITION IS IN GEOGRAPHIC COORDINATES, C AND THE OUTPUT IS IN GEOGRAPHIC COORDINATES. C THE MAGNETIC FIELD AND PLASMA MODELS EXPECT GEOGRAPHIC COORDINATES AS C INPUT. INPUT IS RADIUS IN RJ, THETA(DEG) AND AZIMUTHAL ANGLE (DEG) C C R, THETA, AND PHI ARE THE STANDARD POLAR COORDINATES C INCLUDE 'JUPCOMMONS.FOR' COMMON/BLKMODEL/ENO,ALTI,IP,CEE,PPL,DELB,RPPL,PHIW,ISHEATH CM COMMONS/BLKPARMS/XGSE,YGSE,ZGSE,SINT,COST,SINP,GST,SLONG,SRASN, CM #SDEC,S1,S2,S3,IONS COMMON/BLKTIME/IIYR,IIDAY,DDSECS COMMON/BLKCONSTANTS/BB1,BB2,BB3,AA,BB,ALF1,ALF2,ALF3,GAM1,GAM2, #GAM3,DELTA1,DELTA2,DELTA3,ENSW COMMON/BSMBLK/BIRSM,BITSM,BIPSM,DBIRDR,DBIRDT,DBIRDP,DBITDR, #DBITDT,DBITDP,DBIPDR,DBIPDT,DBIPDP,ITIMES COMMON/BLKOUT/IFRX,WINDOWF COMMON/BLKTOR/CLTORFAC,WMTORFAC,PLASFAC COMMON/BLKDIP/IDIP,THETADIP,PHIDIP,THETMAG,BLD,IMAG CHARACTER*40 FILESTUFF WRITE(6,*) ' INPUT IPLUS: 1=OMODE 2=XMODE' READ(5,*) IPLUS WRITE(6,*) ' 1=O6 2=O4 MAGNETIC FIELD MODEL' READ(5,*) IMAG WRITE(6,*) ' MAXIMUM PLANETARY RADIUS (RJ)' READ(5,*) ALOOK(10) write(6,*) 'INPUT Wave normal angle in degrees' read(5,*) psii NFREQ=1 NPSI=1 PSIINC=0.0 CM PRAD=6378. !EARTH RADIUS IN KM. PRAD=71300. !RADIUS IN KM ALOOKSTP=10 ALOOK(1)=2.0 ALOOK(2)=4.0 ALOOK(3)=6.0 ALOOK(4)=8.0 ALOOK(5)=10.0 ALOOK(6)=12.0 ALOOK(7)=15.0 ALOOK(8)=20.0 ALOOK(9)=25.0 C C INPUT FREQ IN HZ, REL=RADIUS IN RJ TH=THETA(DEG) PHID=PHI(DEG) C WRITE(6,*) ' INPUT F(Hz), SYSIII: R(RJ), THETA(DEG), PHI(DEG)' read(5,*,ERR=10000)FF(1),REL(1),TH(1),PHID(1) WRITE(6,*) ' 1= PLOT FILE; 2=NONE' READ(5,*) IPLT WRITE(6,*) ' 1= .DAT FILE; 2=NONE' READ(5,*) IWRITE WRITE(6,*) ' 1= FRX FILE' READ(5,*) IFRX IF(IFRX.EQ.1) THEN WRITE(6,*) ' INPUT FREQ. WINDOW (KHZ) FOR F-FRX' READ(5,*) WINDOWF ENDIF !IFRX if(th(1).lt.90.0) then pole='N' else pole='S' end if IF(IPLT.EQ.1) THEN OPEN (UNIT=7,NAME='JUPITER2.PLT',TYPE='NEW',CARRIAGECONTROL= #'LIST') ENDIF !IPLT if(iwrite.eq.1) then OPEN (UNIT=8,NAME='JUPITER2.DAT',TYPE='NEW',CARRIAGECONTROL= #'LIST') end if IF(IFRX.EQ.1) THEN OPEN (UNIT=13,NAME='JUPITER2.FRX',TYPE='NEW',CARRIAGECONTROL= #'LIST') WRITE(6,*) ' INPUT COMMENT FOR "JUPITER2.FRX" FILE' READ(5,1339) FILESTUFF WRITE(13,1339) FILESTUFF ENDIF !IFRX 1339 FORMAT(1X,A40) 1304 FORMAT(1X,A70) C RADEG=3.1415927/180. CM WRITE(6,*) 'INPUT IYR,IDAY,DSECS' CM READ(5,*) IYR,IDAY,DSECS C c-------------------------------------------------------------------------- C FOR BR>0 AND BT>0 IN THE NORTHERN HEMISPHERE: C BETA INCREASES IN THE COUNTERCLOCKWISE DIRECTION C AS YOU LOOK INTO B. BETA=90 DEGREES IS IN THE R-THETA PLANE AND C TOWARD THE NORTH POLE. C C FOR BR>0 AND BT>0 IN THE SOUTHERN HEMISPHERE: C BETA=90 DEGREES IS IN THE R-THETA PLANE AND TOWARD THE SOUTH POLE. C WRITE(6,*) ' INPUT NUMBER OF LINES TO SKIP FOR PRINT-OUTS' READ(5,*) IWRIT WRITE(6,*) ' # OF BETAS, INITIAL BETA(DEG), BETA INCREMENT' READ(5,*) NBETA,BETAI,BETINC WRITE(6,*) ' INPUT NUMBER OF IONS' READ(5,*) IONS C WRITE(6,*) 'INPUT COOL TORUS DENSITY FRACTION' READ(5,*) CLTORFAC WRITE(6,*) 'INPUT WARM TORUS DENSITY FRACTION' READ(5,*) WMTORFAC WRITE(6,*) 'INPUT INNER PLASMASPHERE DENSITY FRACTION' READ(5,*) PLASFAC C IF (IFRX.EQ.1) THEN WRITE(13,3373) CLTORFAC,WMTORFAC,PLASFAC 3373 FORMAT(3X,'CL TORUS FAC=',1P1E12.4,2X,'WM TORUS FAC=', # 1P1E12.4,2X,'INPLS FAC=',1P1E12.4) WRITE(13,3377) FF(1),PSII,BETAI 3377 FORMAT(2X,'F=',1P1E12.3,2X,' PSI=',1P1E12.4,3X,'BETA=', #1P1E12.4) ENDIF !IFRX C ICUTLIM=200 C C SET PARAMTERS FOR PRINTOUT IN "OUTP" C pi=-1.0 PI=ACOS(pi) FACTR=1.0 IKEY=0 C "ISTP" IS SET IN ROUTINE FCT WHEN ISTP=1, THEN THE C VALUES OF "Y(4)" AND "Y(5)" ARE TOO LARGE AND THE PROGRAM C IS SENT BACK TO THE MAIN ROUTINE TO START THE NEXT "BETA". ISTP=0 DEGRAD=180.0/PI RADEG=PI/180.0 NEGCNT=0 C C*** NFREQ IS THE NUMBER OF FREQUENCIES, REL'S,TH'S, AND PHI'S. TH IS THETA C NRAD=7 NTYME=KJ*NPSI C C NOW BEGINS THE FREQUENCY LOOP C DO 1001 KK=1,NFREQ F=FF(KK) RJ=REL(KK) THETD=TH(KK) PHIDEG=PHID(KK) THET1=THETD*RADEG PHI1=PHIDEG*RADEG C C*** NPSI IS THE NUMBER OF PSI VALUES C DO JPSI=1,NPSI JJ=1 PSI=PSII+PSIINC*FLOAT(JPSI-1) C SET THE RADIUS PRINT OUT INDEX AND FLAG DO JBETA=1,NBETA ICUT=1 BETA=(BETAI+FLOAT(JBETA-1)*BETINC) IF(BETA.GT.360.) BETA=BETA-360.0 BETAD=BETA BETA=BETA*RADEG IF(IPLT.EQ.1) THEN CM WRITE(7,*)9999.,9999.,9999.,9999.,9999. WRITE(7,1301) F,RJ,THETD,PHIDEG WRITE(7,1302) BETAD,PSI,ITIMES 1301 FORMAT(1X,1P4E12.3) 1302 FORMAT(1X,1P2E12.3,I3) ENDIF !IPLT PSIR=PSI*RADEG pi=-1.0 pi=ACOS(pi) Y1RJ=REL(1) C RAD=57.29578 !DEGREES/RADIAN C CALL MAGNETDRIVER(Y1RJ,THET1,PHI1) C C ALPHA IS THE ANGLE BETWEEN B AND THE PHI DIRECTION C GAMMA IS THE ANGLE BETWEEN THE B FIELD PROJECTED ONTO THE R,THETA PLANE AND R C ALPHA=ACOS(BP/BO) sinalf=SIN(alpha) STEPH=ABS(BR/(BO*SINALF)) GAMMA=ACOS(STEPH) IF((BR.LT.0.0).AND.(BT.GT.0.0)) GAMMA=PI-GAMMA !2ND QUADRANT IF((BR.LT.0.0).AND.(BT.LT.0.0)) GAMMA=GAMMA+PI !3RD QUAD IF((BR.GT.0.0).AND.(BT.LT.0.0)) GAMMA=2.0*PI-GAMMA !4TH QUAD OMEGA=ACOS(COS(ALPHA)*COS(PSIR)+SIN(ALPHA)*SIN(PSIR)* #COS(BETA)) EPSR=pi/2.0 - OMEGA DELNUM=SIN(ALPHA)*SIN(GAMMA)*COS(PSIR)-COS(ALPHA)*SIN(GAMMA)* #SIN(PSIR)*COS(BETA)-COS(GAMMA)*SIN(PSIR)*SIN(BETA) DELDEN=SIN(ALPHA)*COS(GAMMA)*COS(PSIR)-COS(ALPHA)*COS(GAMMA)* #SIN(PSIR)*COS(BETA)+SIN(GAMMA)*SIN(PSIR)*SIN(BETA) RATDEL=ABS(DELNUM/DELDEN) DELR=ATAN(RATDEL) IF((DELDEN.LT.0.0).AND.(DELNUM.GT.0.0)) DELR=PI-DELR !2ND QUAD IF((DELDEN.LT.0.0).AND.(DELNUM.LT.0.0)) DELR=DELR+PI !3RD QUAD IF((DELDEN.GT.0.0).AND.(DELNUM.LT.0.0)) DELR=2.0*PI-DELR !4TH QUAD DO 200 I=1,6 DERY(I)=1.0 200 CONTINUE WRITE(6,213) ALPHA,GAMMA,OMEGA,EPSR 213 FORMAT(1X,'INITIAL ALPHA,GAMMA,OMEGA,EPSR(RAD.) ',1P4E12.3) C C CALCUALTE MODEL CONSTANTS C ZO=0.0 CM IONS=0 W(1)=1.0 !H+ W(2)=16.0 !O+ W(3)=16.0 !O++ W(4)=32.0 !S+ W(5)=32.0 !S++ W(6)=32.0 !S+++ W(7)=23.0 !Na+ ANEO=0.0 C C INITIALIZE PROGRAM (ACC,STEP SIZE,ALTITUDE START&STOP) C ALT1=0.0 ALT2=0.0 ALTUP=(2.0+ALOOK(ALOOKSTP))*PRAD ALTDN=-2.0E-01*PRAD STEP=PRAD/300. STEPMX=PRAD ACC=5.0E-03 IPRINT=1 IC=IPRINT-1 IB=100 IS=0 ITMAX=25 IT=1 ITURN=0 ILINES=0 IBOM=0 C C INITIALIZE PROBLEM C Y(1)=R, Y(2)=COLAT=THETA, Y(3)=PHI, Y(4)=DELR, Y(5)=EPSR, Y(6)=PSIR C C DEFINE THE CHARGE STATE OF IONS USED C ZI(1)=1.0 !ELECTRONS ZI(2)=1.0 !H+ ZI(3)=1.0 !O+ ZI(4)=2.0 !O++ ZI(5)=1.0 !S+ ZI(6)=2.0 !S++ ZI(7)=3.0 !S+++ ZI(8)=1.0 !Na+ ALT=(RJ-1.0)*PRAD Y(1)=PRAD*RJ Y(2)=THETD*RADEG Y(3)=PHIDEG*RADEG Y(4)=DELR Y(5)=EPSR Y(6)=PSIR PRMT(1)=0.0 PRMT(2)=1.2*alook(10)*prad PRMT(3)=STEP PRMT(4)=ACC BETAD=BETA*DEGRAD EPSD=EPSR*DEGRAD DELD=DELR*DEGRAD NDIM=5 JFLG=0 SLIMIT=3.0 CALL DHPCG(PRMT,Y,DERY,NDIM,IHLF,AUX,IB) ISTP=0 END DO !JBETA 102 FORMAT(/,5X,2HF=,F13.2,3X,3HZO=,F7.2,3X,//1X,'ION MASSES= ',7F6.2) 107 FORMAT(1H1) MMMK=JJ-1 IF(MMMK.GT.0) GO TO 120 121 FORMAT(5X,'NO VALUES FOR NRAD PTS.,PSI= ',E12.3) GO TO 999 120 CONTINUE 131 FORMAT(1X,' F(HZ),PSI= ',1P2E12.3) 123 FORMAT(10X,1P4E12.3,2I5) 124 FORMAT(20X,'KKL=',I5,2X,'MMMK=',I5) 130 FORMAT(4(1X,E12.3)) 999 CONTINUE END DO !JPSI 1001 CONTINUE STOP 10000 END C C SUBROUTINE OUTP(X,Y,DERY,IHLF,PRMT,NDIM) C C INCLUDE 'JUPCOMMONS.FOR' COMMON/BLKPROB/OLDR,OLDT,OLDP,OLDY4,OLDY5 COMMON/BLKMODEL/ENO,ALTI,IP,CEE,PPL,DELB,RPPL,PHIW,ISHEATH COMMON/BLKOUT/IFRX,WINDOWF COMMON/BLKDIP/IDIP,THETADIP,PHIDIP,THETMAG,BLD,IMAG CM COMMON/DTEST/OLDDFP2DR,OLDDFP2DT,OLDDFP2DP,DFP2DRS,DFP2DTS, CM #DFP2DPS,DNDXMAG,ANG1,ATANG,DLDXMAG,CTFAC,PP,EL CM COMMON/BSMBLK/BIRSM,BITSM,BIPSM,DBIRDR,DBIRDT,DBIRDP,DBITDR, CM #DBITDT,DBITDP,DBIPDR,DBIPDT,DBIPDP,ITIMES pi=-1.0 pi=ACOS(pi) 999 continue if(y(2).gt.pi) then print *, 'y(2)=',y(2) y(2)=2.0*pi-y(2) goto 999 else if(y(2).lt.0.0) then print *, 'y(2)=',y(2) y(2)=abs(y(2)) goto 999 end if if(y(3).gt.2.0*pi) then print *, 'y(3)=',y(3) y(3)=y(3)-2.0*pi goto 999 else if(y(3).lt.0.0) then print *, 'y(3)=',y(3) y(3)=2.0*pi+y(3) goto 999 end if Y1=Y(1) Y2=Y(2) Y3=Y(3) Y4=Y(4) Y5=Y(5) Y6=Y(6) LAT=90.0-Y2*DEGRAD THETD=Y2*DEGRAD PHIA=Y3*DEGRAD CIRC=360. DELA=MOD(Y4*DEGRAD,CIRC) PSI=Y6*DEGRAD SINY2=SIN(Y2) ALT=Y1-PRAD RJ=Y1/PRAD BL=RJ/(SINY2*SINY2) !L-VALUE ONLY FOR MAGNETIC COORDS. RX=RJ*COS(LAT*RADEG) RY=RJ*SIN(LAT*RADEG) IF(ABS(Y6).GT.10.0) GO TO 300 IF(ALT.GT.ALTUP) GO TO 298 IF(ALT.LT.ALTDN) GO TO 298 GO TO 400 296 CONTINUE 297 FORMAT(10X,'ITURN= ',I5,1X,'ITMAX= ',I5) GO TO 300 298 CONTINUE 299 FORMAT(10X,'ALT= ',E12.4,1X,'ALTUP= ',E12.4,1X,'ALTDN= ',E12.4) 300 CONTINUE PRMT(5)=1.0 GO TO 500 400 CONTINUE IC=IC+1 IF(IC.NE.IPRINT) GO TO 1000 IC=0 500 CONTINUE CM CALL MODEL(Y) FLHR1=0.0 C GF AND PF ARE IN KHZ PF=SQRT(FP2(1))*1.E-3 GF=FB(1)*1.E-3 GFD2=FB(1)/2.0 FRL=SQRT(GFD2*GFD2+FP2(1)) FRO=(GFD2+FRL)*1.E-3 FLO=(-GFD2+FRL)*1.E-3 FKHZ=F*1.E-3 DIFF=ABS(FKHZ-FRO) IF((IFRX.EQ.1).AND.(DIFF.LE.WINDOWF)) THEN WRITE(13,533) RJ,THETD,PHIA,FRO,BLD 533 FORMAT(1X,1P5E12.4) ENDIF !DIFF RES=90.0 IF((-S/P).GT.0.0) RES=RES-DATAN(SQRT(-S/P))*DEGRAD TEV=43.0*BL-29.0 VTH=4.19E07*SQRT(TEV) VP=3.0E10/N ELAMBDA=VP/F CM DEBYE=7.43E02*SQRT(TEV/ANER) CM IF(JFLG.EQ.1) GO TO 601 CM RR=Y(1) CM THAYTA=Y(2) CM PHY=Y(3) CM DIST=0.0 CM STIME=0.0 CM601 DR=Y(1)-RR CM DTHATA=Y(2)-THAYTA CM DPHY=Y(3)-PHY CM625 FORMAT(2X,'ENGR=1',2X,'PF=',E12.4,2X,'GF=',E12.4) IF(MOD(ILINES,100).NE.0) GO TO 701 IF(IWRITE.EQ.1) WRITE(8,339) BETAD 339 FORMAT(2X,'BETA=',1P1E12.3) IF(IWRITE.EQ.1) WRITE(8,700) IF(IWRITE.EQ.1) WRITE(8,733) IBOM=IBOM+1 IF(IBOM.GT.500) PRMT(5)=1.0 702 FORMAT(10X,'IBOM IS > 500') 700 FORMAT(6X,'RJ',11X,'BL',8X,'THETA',8X,'PHI',9X,'PSI',9X,'ANER', #7X,'FP(KHZ)',5X,'FG(KHZ)',6X,' N ',' VGROUP ',' ALPHA') 733 FORMAT(1X,' FRX ',' FLO ',' BR ', #' BT ',' BP ',' GAMMA ',' ALPHA ', #' OMEGA ',' DEL ',' POLAR ') 701 CONTINUE YPR2D=YPR2*DEGRAD C C NEW SCEME FOR PRINTOUT C XGSEE=RJ*SIN(Y2)*COS(Y3) YGSEE=RJ*SIN(Y2)*SIN(Y3) ZGSEE=RJ*COS(Y2) SGSEE=SQRT(XGSEE*XGSEE+YGSEE*YGSEE+ZGSEE*ZGSEE) DXGSE=OLDXGSEE-XGSEE DYGSE=OLDYGSEE-YGSEE DZGSE=OLDZGSEE-ZGSEE DSGSE=SQRT(DXGSE*DXGSE+DYGSE*DYGSE+DZGSE*DZGSE) DSGSEKM=DSGSE*PRAD IF(VGROUP.NE.0.0) THEN VGROUPKM=VGROUP*3.0E5 DTGSEKM=DSGSEKM/VGROUPKM DELAYT=DELAYT+DTGSEKM !TOTAL DELAY TIME (SEC) ELSE !VGROUP=0 DELAYT=1000.0 ENDIF !VGROUP OLDXGSEE=XGSEE OLDYGSEE=YGSEE OLDZGSEE=ZGSEE IF(MOD(ILINES,IWRIT).NE.0) GO TO 713 C if(iplt.eq.1) WRITE(7,1303) XGSEE,YGSEE,ZGSEE,SGSEE,POLAR, # DELAYT,VGROUP CM if(iplt.eq.1)WRITE(7,1303) RJ,THETD,PHIA,DELA,BETAD 1303 FORMAT(1X,1P7E12.4) IF(IWRITE.EQ.1) THEN WRITE(8,800) RJ,EL,THETD,PHIA,PSI,ANER,PF,GF,CN,VGROUP, #ALFSTIX WRITE(8,803) FRO,FLO,BR,BT,BP,GAMMAD,ALPHAD,OMEGAD,DELD, #POLAR CM WRITE(8,803) OLDDFP2DR,DFP2DR(1),OLDDFP2DT,DFP2DT(1),OLDDFP2DP, CM #DFP2DP(1),DFP2DRS,DFP2DTS,DFP2DPS CM WRITE(8,804) DNDXMAG,ATANG,DLDXMAG,ANG1,PP CM WRITE(8,800) RE,BL,THETD,PHIA,PSI,ANER,PF,GF,N,DIST,STIME,VGR, CM #VP,DNEDR END IF !IWRITE C######### C C SCHEME FOR PREVENTING A RAY FROM BEING STUCK AT ONE POSITION C C######### 713 ILINES=ILINES+1 CM WRITE(6,*) ' ICUT,ICUTLIM',ICUT,ICUTLIM OLDR=Y(1) OLDT=Y(2) OLDP=Y(3) OLDY4=Y(4) OLDY5=Y(5) DIFRJ=ABS(OLDRJ-RJ) CM DIFR=ABS(OLDR-Y(1)) IF(DIFRJ.LT.1.0D-06) THEN ICUT=ICUT+1 IF(ICUT.GE.ICUTLIM) THEN if(iwrite.eq.1) then WRITE(8,*) ' VECTOR IS STUCK; INCREMENT BETA' WRITE(8,*) ' RJ,OLDRJ,ICUT= ',RJ,OLDRJ,ICUT end if PRMT(5)=1.0 END IF !ICUT > ICUTLIM ELSE !DIFRJ ICUT=1 ENDIF !DIFRJ OLDRJ=RJ 800 FORMAT(1X,1P11E12.4) 803 FORMAT(2X,1P10E12.4) 804 FORMAT(2X,1P5E12.4) CM800 FORMAT(1X,1P11E12.4,/,2X,1P3E12.4) IF (PRMT(5).NE.1.0) GO TO 1000 900 FORMAT(/,2X,3HRJ=,F6.2,4X,3HRX=,F6.2,4X,3HRY=,F6.2,3X,4HLAT=,F6.2 1 ,4X,5HPHID=,F7.2,/,2X,8HFP(KHZ)=,F8.2,3X,8HFG(KHZ)=,F8.2,3X,4HPSI= 2 ,F7.2,3X,5HCCHI=,F7.2,///) 1000 CONTINUE JFLG=1 RETURN END C FUNCTION XTAN2(Y,X) IMPLICIT real*8 (A-H,O-Z) A = ABS(X) B = ABS(Y) IF(X.EQ.0.0) GO TO 40 C = B/A C = DATAN(C) IF(X.LT.0.0) GO TO 10 IF(Y.LT.0.0) GO TO 30 XTAN2 = C RETURN 10 CONTINUE IF(Y.LT.0.0) GO TO 20 xtan2=-1.0 XTAN2 = ACOS(xtan2) -C RETURN 20 CONTINUE xtan2=-1.0 XTAN2 = C - ACOS(xtan2) RETURN 30 CONTINUE XTAN2=-C RETURN 40 CONTINUE xtan2=1.0 XTAN2 = ASIN(xtan2) IF(Y.LT.0.0) XTAN2 =-XTAN2 RETURN END C C SUBROUTINE DHPCG(PRMT,Y,DERY,MDIM,IHLF,AUX,IB) C MOFIFIED 09/19/70 TO UPDATE VALUES OF Y FROM PREDICTOR AND CORRECTOR C MODIFIED 09/21/72 TO CHECK ACCURACY OF ALL PARAMETERS C MODIFIED 2/1/73 TP CHECK WEIGHTED PERCENTANGE ACCURACY OF ALL PARAMERERS C INCLUDE 'JUPCOMMONS.FOR' COMMON/BLKPROB/OLDR,OLDT,OLDP,OLDY4,OLDY5 NDIM=MDIM IHLF=0 C C COUNTER FOR ISTP FLAG C xlim=prad/2.0 CM hlim=prad/100.0 hlim=prad/20.0 XMXLIM=10.*PRAD ISTPCNT=1 N=1 X=PRMT(1) H=PRMT(3) 1313 IF(ISTP.EQ.1) THEN H=H/2.0 Y(1)=OLDR Y(2)=OLDT Y(3)=OLDP Y(4)=OLDY4 Y(5)=OLDY5 IF(IWRITE.EQ.1) WRITE(8,*) ' D-IF: Y1= ',Y(1) ISTP=0 ISTPCNT=ISTPCNT+1 IF(ISTPCNT.GT.20) THEN if(iwrite.eq.1) WRITE(8,*) ' ISTPCNT > 20' ISTPCNT=1 RETURN ENDIF !ISTPCNT ENDIF !ISTP PRMT(5)=0.0 DO 1 I=1,NDIM AUX(16,I)=0.0 AUX(15,I)=DERY(I) AUX(1,I)=Y(I) 1 CONTINUE IF(H*(PRMT(2)-X))3,2,4 2 IHLF=50 GOTO 4 3 IHLF=60 4 CALL FCT(X,Y,DERY,2) IF(ISTP.EQ.1) THEN GO TO 1313 ENDIF !ISTP CALL OUTP(X,Y,DERY,IHLF,PRMT,NDIM) IF(PRMT(5))6,7,6 6 WRITE(6,*) 'PRMT=1 LABELED LINE 6' RETURN 7 DO 8 I=1,NDIM AUX(8,I)=DERY(I) 8 CONTINUE ISW=1 GOTO 100 9 X=X+H CALL FCT(X,Y,DERY,1) IF(ISTP.EQ.1) GO TO 1313 DO 10 I=1,NDIM AUX(2,I)=Y(I) 10 CONTINUE 11 IHLF=IHLF+1 X=X-H DO 12 I=1,NDIM AUX(4,I)=AUX(2,I) 12 CONTINUE H=H/2.0 if(x.lt.xlim) then CM if((x.lt.xlim).OR.(X.GT.XMXLIM)) then if(h.gt.hlim) h=hlim end if N=1 ISW=2 GOTO 100 13 X=X+H CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 N=2 DO 14 I=1,NDIM AUX(2,I)=Y(I) AUX(9,I)=DERY(I) 14 CONTINUE ISW=3 GOTO 100 15 DELT=0.0 IF(IFLAG.EQ.1) GO TO 17 DO 16 I=1,NDIM IF( ABS(Y(I)-AUX(4,I))*AUX(15,I).GT.PRMT(4)* ABS(Y(I))*15.0) 1 GO TO 17 16 CONTINUE GO TO 19 17 IF(IHLF-IB)11,18,18 18 IHLF=99 X=X+H GOTO 4 19 X=X+H CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 DO 20 I=1,NDIM AUX(3,I)=Y(I) AUX(10,I)=DERY(I) 20 CONTINUE N=3 ISW=4 GOTO 100 21 N=1 X=X+H CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 X=PRMT(1) DO 22 I=1,NDIM AUX(11,I)=DERY(I) Y(I)=AUX(1,I)+H*(.375*AUX(8,I)+.79166667*AUX(9,I) - 1 0.20833334*AUX(10,I)+.04166667*DERY(I)) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D1: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 22 CONTINUE 23 X=X+H N=N+1 CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 CALL OUTP(X,Y,DERY,IHLF,PRMT,NDIM) IF(PRMT(5))6,24,6 24 IF(N-4)25,200,200 25 DO 26 I=1,NDIM AUX(N,I)=Y(I) AUX(N+7,I)=DERY(I) 26 CONTINUE IF(N-3)27,29,200 27 DO 28 I=1,NDIM DELT=AUX(9,I)+AUX(9,I) DELT=DELT+DELT Y(I)=AUX(1,I)+.33333333*H*(AUX(8,I)+DELT+AUX(10,I)) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D2: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 28 CONTINUE GOTO 23 29 DO 30 I=1,NDIM DELT=AUX(9,I)+AUX(10,I) DELT=DELT+DELT+DELT Y(I)=AUX(1,I)+.375*H*(AUX(8,I)+DELT+AUX(11,I)) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D3: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 30 CONTINUE GOTO 23 100 DO 101 I=1,NDIM Z=H*AUX(N+7,I) AUX(5,I)=Z Y(I)=AUX(N,I)+.4*Z IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D4: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 101 CONTINUE Z=X+.4*H CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 DO 102 I=1,NDIM Z=H*DERY(I) AUX(6,I)=Z Y(I)=AUX(N,I)+.29697761*AUX(5,I)+.15875964*Z IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D5: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 102 CONTINUE Z=X+.45573725*H DO 103 I=1,NDIM Z=H*DERY(I) AUX(7,I)=Z Y(I)=AUX(N,I)+.21810039*AUX(5,I)-3.0509651*AUX(6,I)+3.8328648*Z IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D6: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 103 CONTINUE Z=X+H CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 DO 104 I=1,NDIM Y(I)=AUX(N,I)+.17476028*AUX(5,I)-.55148066*AUX(6,I)+ 1 1.2055356*AUX(7,I)+.17118478*H*DERY(I) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D7: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 104 CONTINUE GOTO(9,13,15,21),ISW 200 ISTEP=3 201 IF(N-8)204,202,204 202 DO 203 N=2,7 DO 203 I=1,NDIM AUX(N-1,I)=AUX(N,I) AUX(N+6,I)=AUX(N+7,I) 203 CONTINUE N=7 204 N=N+1 DO 205 I=1,NDIM AUX(N-1,I)=Y(I) AUX(N+6,I)=DERY(I) 205 CONTINUE X=X+H 206 ISTEP=ISTEP+1 DO 207 I=1,NDIM DELT=AUX(N-4,I)+1.3333333*H*(AUX(N+6,I)+AUX(N+6,I)- 1 AUX(N+5,I)+AUX(N+4,I)+AUX(N+4,I)) CM207 Y(I)=DELT Y(I)=DELT IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D8: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 207 CONTINUE CALL FCT(X,Y,DERY,1) IF(ISTP.EQ.1) GO TO 1313 DO 2070 I=1,NDIM DELT=Y(I) Y(I)=DELT-.92561983*AUX(16,I) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D9: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 AUX(16,I)=DELT 2070 CONTINUE C IF(Y(1).LT.0.0) THEN ISTP=1 GO TO 1313 ENDIF !Y1 < 0 CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 DO 208 I=1,NDIM DELT=.125*(9.0*AUX(N-1,I)-AUX(N-3,I)+3.0*H*(DERY(I)+AUX(N+6,I)+ 1 AUX(N+6,I)-AUX(N+5,I))) CM208 Y(I)=DELT Y(I)=DELT IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D10: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 208 CONTINUE CALL FCT(X,Y,DERY,1) IF(ISTP.EQ.1) GO TO 1313 DO 2080 I=1,NDIM DELT=Y(I) AUX(16,I)=AUX(16,I)-DELT Y(I)=DELT+.07438017*AUX(16,I) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D11: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 2080 CONTINUE DELT=0.0 IF(IFLAG.EQ.1) GO TO 222 DO 209 I=1,NDIM IF( ABS(AUX(16,I))*AUX(15,I).GT.PRMT(4)* ABS(Y(I))) GO TO 222 209 CONTINUE 210 CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 CALL OUTP(X,Y,DERY,IHLF,PRMT,NDIM) IF(PRMT(5))212,211,212 211 IF(IHLF-IB-1) 213,212,212 212 if(iwrite.eq.1)WRITE(8,2121) X,H,PRMT(5),IHLF,IB 2121 FORMAT(2X,'X=',2X,E12.4,2X,'H=',2X,E12.4,2X,'PRMT5=',2X,E12.4, 1 'IHLF=',2X,I4,2X,'IB=',2X,I4,/1X,'DHPCG,212') RETURN 213 IF(H*(X-PRMT(2)))214,555,555 555 if(iwrite.eq.1)WRITE(8,556) H,X,PRMT(2) 556 FORMAT('H= ',E12.3,1X,'X= ',E12.3,1X,'PRMT2= ',E12.3,2X,'DHPCG,213') RETURN 214 IF( ABS(X-PRMT(2))-.1D0* ABS(H))555,215,215 215 DO 227 I=1,NDIM IF(1.E1* ABS(AUX(16,I))*AUX(15,I).GT.PRMT(4)* ABS(Y(I)))GOTO201 227 CONTINUE 217 IF(N-7)201,218,218 218 IF(ISTEP-4)201,219,219 219 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)201,220,201 220 H=H+H IF(H.GT.STEPMX) H=STEPMX if(x.lt.xlim) then CM if((x.lt.xlim).OR.(X.GT.XMXLIM)) then if(h.gt.hlim) h=hlim end if IHLF=IHLF-1 ISTEP=0 DO 221 I=1,NDIM AUX(N-1,I)=AUX(N-2,I) AUX(N-2,I)=AUX(N-4,I) AUX(N-3,I)=AUX(N-6,I) AUX(N+6,I)=AUX(N+5,I) AUX(N+5,I)=AUX(N+3,I) AUX(N+4,I)=AUX(N+1,I) DELT=AUX(N+6,I)+AUX(N+5,I) DELT=DELT+DELT+DELT AUX(16,I)=8.9629630*(Y(I)-AUX(N-3,I))-3.3611111*H*(DERY(I)+ 1 DELT+AUX(N+4,I)) 221 CONTINUE GOTO 201 222 IHLF=IHLF+1 IF(IHLF-IB)223,223,210 223 H=H/2.0 if(x.lt.xlim) then CM if((x.lt.xlim).OR.(X.GT.XMXLIM)) then if(h.gt.hlim) h=hlim end if ISTEP=0 DO 224 I=1,NDIM Y(I)=.390625E-2*(8.E1*AUX(N-1,I)+135.0*AUX(N-2,I)+4.E1*AUX(N-3,I) 1 +AUX(N-4,I))-.1171875*(AUX(N+6,I)-6.0*AUX(N+5,I)-AUX(N+4,I))*H AUX(N-4,I)=.390625E-2*(12.0*AUX(N-1,I)+135.0*AUX(N-2,I)+ 1 108.0*AUX(N-3,I)+AUX(N-4,I))-.0234375*(AUX(N+6,I)+ 2 18.0*AUX(N+5,I)-9.0*AUX(N+4,I))*H AUX(N-3,I)=AUX(N-2,I) AUX(N+4,I)=AUX(N+5,I) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D12: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 224 CONTINUE X=X-H DELT=X-(H+H) CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 DO 225 I=1,NDIM AUX(N-2,I)=Y(I) AUX(N+5,I)=DERY(I) Y(I)=AUX(N-4,I) IF(Y(1).LT.0.0) THEN IF(IWRITE.EQ.1) WRITE(8,*) ' D13: Y1,H,P2= ',Y(1),H,PRMT(2) ENDIF !Y1 < 0 225 CONTINUE DELT=DELT-(H+H) CALL FCT(X,Y,DERY,0) IF(ISTP.EQ.1) GO TO 1313 DO 226 I=1,NDIM DELT=AUX(N+5,I)+AUX(N+4,I) DELT=DELT+DELT+DELT AUX(16,I)=8.9629630*(AUX(N-1,I)-Y(I))- 1 3.3611111*H*(AUX(N+6,I)+DELT+DERY(I)) AUX(N+3,I)=DERY(I) 226 CONTINUE GOTO 206 END C SUBROUTINE FCT (X,Y,DERY,NUP) C INCLUDE 'JUPCOMMONS.FOR' CM IMPLICIT real*8 (A-H,O-Z) CM real*8 L,NR,NT,NP CM REAL*8 LAT,LONG CM character*1 pole,kind CM DIMENSION FF(10),REL(10),TH(10),PHID(10),ALOOK(10),AUX(16,6) CM DIMENSION Y(6),DERY(6),PRMT(5),ZI(8),W(8) CM DIMENSION FP2(8),DFP2DR(8),DFP2DT(8),DFP2DP(8),PLASCON(8) CM DIMENSION FB(8),FB2(8),DFBDR(8),DFBDT(8),DFBDP(8) CM DIMENSION GG(8),E(8),DRPDFP(8),DLDFP(8), CM 2 DLPDFP(8),DPPDFP(8),DLPDFB(8),DPPDFB(8), CM 3 DPDFP(8),DRPDFB(8),DLDFB(8),DAPDFP(8),DBPDFP(8),DCPDFP(8), CM 4 DN2DFP(8),DAPDFB(8),DBPDFB(8),DCPDFB(8),DN2DFB(8) CM DIMENSION DNDFP(8),DNDFB(8) C CM COMMON /INDEX/CN,S,P,D,SIP,COP,ALTUP,ALTDN,F,ALT,ACC, CM 1 COSDEL,SINDEL,RO,W,NR,NT,NP,FP,FG,T,ALT1, CM 2 ALT2,IC,IONS,IPRINT,ITURN,ITMAX,ILINES,IT CM common/index2/RE,PRAD,JK,ISERVO,STEP CM COMMON /BLK001/ DEGRAD,RADEG,BR,BT,BP,BO,DBODR,DBODT,DBODP, CM 1 DBRDR,DBRDT,DBRDP,DBTDR,DBTDT,DBTDP,STEPMX CM COMMON /JUNK4/BETAD,CHKB,CHKR,CHKL,CHKP,JFLG,NRAD,ALOOK, CM 1 ALOOKSTP,SLIMIT CM COMMON /JUNK5/FACTR,jj,ibom,iflag CM COMMON /MAIN/NEGCNT,IKEY,IWRITE CM COMMON/FFCT/ISTP CM COMMON/BLKA1/ZI CM COMMON/BLKCUT/OLDRJ,ICUTLIM,ICUT,IWRIT,iplt CM COMMON/BLKSIGN/IPLUS,ISIGN,mflag,nflag CM COMMON /JUNK2/YPR2,pole,kind CM COMMON/GYRO/FB,FB2,DFBDR,DFBDT,DFBDP CM COMMON/PLAS/FP2,DFP2DR,DFP2DT,DFP2DP,ANER CM COMMON/DGPLAS/PLASCON COMMON/BLKPROB/OLDR,OLDT,OLDP,OLDY4,OLDY5 real*8 L1,N,M2,M4,LP,L1PRM DATA E(1)/-1.0/,E(2)/1.0/,E(3)/1.0/,E(4)/1.0/ IFLAG=0 C C** THESE LINES PREVENT PHI FROM BEING >2PI OR <0. C twopi=-1.0 pi=ACOS(twopi) TWOPI=2.0*pi Y(3)=MOD(Y(3),TWOPI) IF(Y(3).LT.0.) Y(3)=TWOPI+Y(3) IF(Y(1).LT.0.0) THEN WRITE(8,*) ' F1: Y(1)= ',Y(1) ENDIF !Y1 < 0 CALL MODEL(Y) IF(ISTP.EQ.1.0) RETURN IONS1=IONS+1 DELR=Y(4) EPSR=Y(5) IF(ABS(Y(4)).LT.25.0) GO TO 49 IFLAG=1 !* IF(IWRITE.EQ.1) WRITE(8,*) 'IFLAG=1; Y(4)= ',Y(4) RETURN 49 IF(ABS(Y(5)).LT.25.0) GO TO 48 IFLAG=1 !* IF(IWRITE.EQ.1) WRITE(8,*) 'IFLAG=1; Y(5)= ',Y(5) RETURN 48 CONTINUE DELD=DELR*DEGRAD OMEGA=90.0*RADEG-EPSR C COSDEL IS THE PROJECTION OF K ONTO THE RADIUS VECTOR R C SINDEL IS THE PROJECION OF K IN THE COLAT DIRECTION COSDEL=COS(DELR) SINDEL=SIN(DELR) COSEPS=COS(EPSR) SINEPS=SIN(EPSR) COSO=COS(OMEGA) SINO=SIN(OMEGA) C ALPHA and GAMMA are magnetic field orientation angles in frame of ray COSA=BP/BO SINA=SQRT(BR**2+BT**2)/BO COSG=BR/(BO*SINA) SING=BT/(BO*SINA) IF(COSG.GT.1.0) COSG=1.0 IF(COSG.LT.-1.0) COSG=-1.0 IF(SING.GT.1.0) SING=1.0 IF(SING.LT.-1.0) SING=-1.0 IF(COSA.GT.1.0) COSA=1.0 IF(COSA.LT.-1.0) COSA=-1.0 if(abs(cosa).lt.1.0d-14) then if(cosa.lt.0.0) then cosa=-1.0d-14 else cosa=1.0d-14 end if end if IF(SINA.GT.1.0) SINA=1.0 IF(SINA.LT.-1.0) SINA=-1.0 CM IF(SING.GE.0.0) THEN CM GAMMA=ACOS(COSG) CM ELSE CM GAMMA=2.0*PI-ACOS(COSG) CM ENDIF ALPHA=ACOS(BP/BO) SINALF=SIN(ALPHA) STEPH=ABS(BR/(BO*SINALF)) IF(STEPH.LT.1) THEN GAMMA=ACOS(STEPH) ELSE GAMMA=0.0 WRITE(6,3322) BR,BT,BP,BO 3322 FORMAT(2X,'FCT: BR,BT,BP,BO',1P4E12.3) ENDIF !STEPH < 1 IF((BR.LT.0.0).AND.(BT.GT.0.0)) GAMMA=PI-GAMMA !2ND QUADRANT IF((BR.LT.0.0).AND.(BT.LT.0.0)) GAMMA=GAMMA+PI !3RD QUAD IF((BR.GT.0.0).AND.(BT.LT.0.0)) GAMMA=2.0*PI-GAMMA !4TH QUAD GAMMAD=GAMMA*DEGRAD ALPHAD=ALPHA*DEGRAD OMEGAD=OMEGA*DEGRAD PSIR=SINO*COSDEL*SINA*COSG+SINO*SINDEL*SINA*SING+COSO*COSA if(psir.gt.1.0) then psir=0.0 else if(psir.lt.-1.0) then psir=pi else PSIR=ACOS(psir) end if Y(6)=PSIR SIP=SIN(PSIR) if(abs(sip).lt.1.0d-10) sip=1.0d-10 COP=COS(PSIR) SIN2=SIP*SIP COS2=COP*COP C C SPECIES LOOP C WF=TWOPI*F WF2=WF*WF DRDWF1=0.0 DRDWF2=0.0 DPDWF1=0.0 C C DEFINE VALUES FOR VGROUP (SHAWHAN, P. 124, 125) C C R,L,P AND REFRACTIVE INDEX C F2=F*F PLASMA=1.0 L1=1.0 R1=1.0 RT1=1.0 DO J=1,IONS1 C ATOP=2.*WF+E(J)*FB(J) ABOT1=(WF+E(J)*TWOPI*FB(J))**2 ABOT2=(WF+E(J)-TWOPI*FB(J))**2 FRATIO=FP2(J)/F2 DRDWF=DRDWF1+FRATIO*ATOP/ABOT1 DLDWF=DLDWF1+FRATIO*ATOP/ABOT2 DPDWF=DPDWF1+2.0*FRATIO/WF C R1=R1-FP2(J)/(F*(F+E(J)*FB(J))) !OLD R VALUE L1=L1-FP2(J)/(F*(F-E(J)*FB(J))) !OLD L VALUE PLASMA=PLASMA-FP2(J)/(F2) !OLD P VALUE END DO C Rs=R1 L=L1 P=P1 D=0.5*(Rs-L) S=0.5*(Rs+L) SAVE=S P=PLASMA CM WRITE(6,53) R,L,S,D CM 53 FORMAT(1X,'R,L,S,D',1P4E12.3) CM WRITE(6,55) FP2(1),FB(1) CM 55 FORMAT(2X,'FP2,FB',1P2E12.3) C A=S*SIN2+P*COS2 B=RS*L*SIN2+P*S*(1.0+COS2) C=P*RS*L DADWF=0.5*(DRDWF+DLDWF)*SIN2+DPDWF*COS2 DBDWF=(RS*DLDWF+L*DRDWF)*SIN2+(0.5*P*(DRDWF+DLDWF)+S*DPDWF)* # (1.0+COS2) DCDWF=P*(RS*DLDWF+L*DRDWF)+RS*L*DPDWF C CCC R,L,P AND REFRACTIVE INDEX C PLASMA=1.0 L1=1.0 C C INSERT NEW PRIME VALUES C YY2=(FB(1)/F)**2 YY1=1.0-YY2 R1PRM=YY1 PLSPRM=YY1 L1PRM=YY1 C DO 100 J=1,IONS1 IF(J.EQ.1) THEN R1PRM=R1PRM-(FP2(1)*(1.0+FB(1)/F))/F2 !NEW PRIMED R L1PRM=L1PRM-(FP2(1)*(1.0-FB(1)/F))/F2 !NEW PRIMED L PLSPRM=PLSPRM-FP2(1)*YY1/F2 !NEW P ELSE CM IF(IWRITE.EQ.1) WRITE(8,*) ' TEST OF FCT J>1 ' R1PRM=R1PRM-(FP2(J)*YY1)/(F*(F+E(J)*FB(J))) !NEW R L1PRM=L1PRM-(FP2(J)*YY1)/(F*(F-E(J)*FB(J))) !NEW L PLSPRM=PLSPRM-FP2(J)*YY1/F2 !NEW P C R1=R1-FP2(J)/(F*(F+E(J)*FB(J))) !OLD R VALUE END IF !J=1 L1=L1-FP2(J)/(F*(F-E(J)*FB(J))) !OLD L VALUE PLASMA=PLASMA-FP2(J)/(F2) !OLD P VALUE 100 CONTINUE RP=R1PRM L=L1 LP=L1PRM DP=0.5*(RP-LP) SP=0.5*(RP+LP) P=PLASMA PPRM=PLSPRM CM WRITE(6,58) RP,LP,DP CM 58 FORMAT(1X,'RP,LP,DP',1P3D12.3) C C DEFINE PRIMED VALUES C APRM=SP*SIN2+PPRM*COS2 BPRM=RP*L*SIN2+P*SP*(1.0+COS2) CPRM=P*L*RP fmult=rp*l-sp*p if(abs(fmult).gt.1.0d14) then ffreq=abs(fmult*sip*sip) else FFreq=SQRT((fmult)**2*SIP**4+4.0*(P*DP)**2*COS2) end if IF(IPLUS.EQ.1) THEN CN2=(BPRM+FFreq)/(2.0*APRM) ELSE CN2=(BPRM-FFreq)/(2.0*APRM) IF(NRITE.EQ.1) THEN RATIOF=F/FB(1) END IF !NRITE END IF !IPLUS IF(CN2.GT.0.0) GO TO 300 CN2=1.0 CN=1.0 IF(IWRITE.EQ.1) WRITE(8,200) 200 FORMAT(2X,'ANOTHER BISECTION WAS MADE TO AVOID IMAGINARY N') ILINES=ILINES+1 IFLAG=1 RETURN 300 CONTINUE C C EVALUATE POLARIZATION C CM DT=DP/YY1 CM if(dt.eq.0.0) dt=1.0d-20 CM ST=SP/YY1 CM POLAR=(CN2-ST)/DT IF(D.EQ.0.0) THEN D=1.E-20 ENDIF !D=0 POLAR=(CN2-S)/D CN=SQRT(CN2) N=CN M2=CN2 M4=M2*M2 C C DEFINE NEW PRIMED DERIVATIVES C CCC DERIVATIVES OF R,L,P WRT FP2,FB,V DO 400 J=1,IONS1 DLDFP(J)=-1.0/(F*(F-E(J)*FB(J))) DPDFP(J)=-1.0/F2 DLPDFP(J)=-YY1/(F*(F-E(J)*FB(J))) DPPDFP(J)=-YY1/F2 IF(J.EQ.1) THEN DRPDFP(1)=-(1.0+FB(1)/F)/F2 DRPDFB(1)=-2.0*FB(1)/F2-FP2(1)/(F2*F) DLPDFB(1)=-2.0*FB(1)/F2+FP2(1)/(F2*F) DPPDFB(1)=-2.0*FB(1)*P/F2 ELSE DRPDFP(J)=-YY1/(F*(F+E(J)*FB(J))) DRPDFB(J)=YY1*FP2(J)*E(J)/(F*(F+E(J)*FB(J))**2) DLPDFB(J)=-YY1*FP2(J)*E(J)/(F*(F-E(J)*FB(J))**2) DPPDFB(J)=0.0 END IF ! J=1 frat=fp2(j)/f DLDFB(J)=E(J)*FRAT/((F+E(J)*FB(J))**2) 400 CONTINUE CON1=-(2.0*APRM*M2-BPRM) IF(CON1.NE.0.) GO TO 401 CON1=1.0D-10 401 CONTINUE CON2=(2.0*M2*(SP-PPRM)*COP-2.0*(RP*L-P*SP)*COP)/CON1 CCC DERIVATIVE OF N WRT FP2,FB,V C DO 500 J=1,IONS1 DAPDFP(J)=.5*(DRPDFP(J)+DLPDFP(J))*SIN2+DPPDFP(J)*COS2 DBPDFP(J)=(RP*DLDFP(J)+L*DRPDFP(J))*SIN2+(P*0.5*(DRPDFP(J)+ 1 DLPDFP(J))+SP*DPDFP(J))*(1.0+COS2) DCPDFP(J)=P*(RP*DLDFP(J)+L*DRPDFP(J))+L*RP*DPDFP(J) DN2DFP(J)=(DAPDFP(J)*M4-DBPDFP(J)*M2+DCPDFP(J))/CON1 DAPDFB(J)=.5*(DRPDFB(J)+DLPDFB(J))*SIN2+DPPDFB(J)*COS2 DBPDFB(J)=(RP*DLDFB(J)+L*DRPDFB(J))*SIN2+P*0.5*(DRPDFB(J)+ 1 DLPDFB(J))*(1.0+COS2) DCPDFB(J)=P*(RP*DLDFB(J)+L*DRPDFB(J)) DN2DFB(J)=(DAPDFB(J)*M4-DBPDFB(J)*M2+DCPDFB(J))/CON1 500 CONTINUE DN2DFB1=DN2DFB(1) DN2DFP1=DN2DFP(1) DFBDR1=DFBDR(1) DFBDT1=DFBDT(1) DFBDP1=DFBDP(1) DFP2DR1=DFP2DR(1) CCC SPATIAL DERIVATIVES OF N Y1=Y(1) if(abs(y1).lt.1.0d-10) y1=1.0d-10 SINT=SIN(Y(2)) if(abs(sint).lt.1.0d-10) sint=1.0d-10 COST=COS(Y(2)) DN2DPS=M2*SIP*CON2 C DN2DX = 2*N*DNDX DNDPS=DN2DPS/(2.0*CN) DNDPS1=DNDPS DN2DPS1=DN2DPS TANALF=-DNDPS/N !THIS IS TAN(ALPHA) FROM STIX (EQN 3-19; P. 53) NR=N*SINO*COSDEL NT=N*SINO*SINDEL NP=N*COSO if(nr.ge.0.0) then if(nr.lt.1.0d-6) nr=1.0d-6 else if(abs(nr).lt.1.0d-6) nr=-1.0d-6 end if if(nt.ge.0.0) then if(nt.lt.1.0d-6) nt=1.0d-6 else if(abs(nt).lt.1.0d-6) nt=-1.0d-6 end if if(np.ge.0.0) then if(np.lt.1.0d-6) np=1.0d-6 else if(abs(np).lt.1.0d-6) np=-1.0d-6 end if IF(EPSR.EQ.0.0) NP=0.0 IF(N.EQ.0) N=1.0D-06 IF(SIP.EQ.0) SIP=1.0D-10 C Spatial derivatives (R,T,P) of PSI needed for calculation of C DNR(,T,P)/DS or DERF,DERG,DERH DPSIDG=SIN(GAMMA-DELR)*COSEPS*SINA/SIP DPSIDA=(-COS(GAMMA-DELR)*COSEPS*COSA+SINEPS*SINA)/SIP DGDR=((COSG/SINA)*DBTDR-(SING/SINA)*DBRDR)/BO DGDT=((COSG/SINA)*DBTDT-(SING/SINA)*DBRDT)/BO DGDP=((COSG/SINA)*DBTDP-(SING/SINA)*DBRDP)/BO DADR=((SING/COSA)*DBTDR+(COSG/COSA)*DBRDR-(SINA/COSA)*DBODR)/BO DADT=((SING/COSA)*DBTDT+(COSG/COSA)*DBRDT-(SINA/COSA)*DBODT)/BO DADP=((SING/COSA)*DBTDP+(COSG/COSA)*DBRDP-(SINA/COSA)*DBODP)/BO DPSIDR=DPSIDG*DGDR+DPSIDA*DADR DPSIDT=DPSIDG*DGDT+DPSIDA*DADT DPSIDP=DPSIDG*DGDP+DPSIDA*DADP 314 FORMAT(1X,'NR',1P4E12.4,/,1P4E12.4) 313 FORMAT(1X,'DPSI',1P6E12.4) C C RAY = N/SQRT(N**2+DNDPS**2) C RAY = 1/(V GROUP*N) if(abs(dn2dps).ge.1.0d+14) then if(dn2dps.gt.0.0) then dn2dps=1.0d+14 else dn2dps=-1.0d+14 end if end if RAY=2.0*M2/SQRT(4.0*m4+dn2dps**2) C VRAY=1.0/(CN*RAY) C C DEFINE VGROUP VIA BUDDEN, PP. 129-131 C C FIRST DNDWF ACCORDING TO SHAWHAN, p. 124, EQN. A40 C DENOMCHK=4.*A*CN*CN2-2.*B*CN IF(DENOMCHK.NE.0.0) THEN DNDWF=-(CN2*CN2*DADWF-CN2*DBDWF+DCDWF)/ # DENOMCHK C ALFSTIX=ABS(ATAN(TANALF)) CNGRP=CN+WF*DNDWF !BUDDEN, EQN. 5.75, P. 131 VGROUP=1.0/(CNGRP*COS(ALFSTIX)) ELSE !(DENOMCHK=0) VGROUP=0.0 ENDIF !DENOMCHK C C DS/ST = VRAY (SEE SHAWHAN, p. 21 AND p. 27) C S = THE RAY PATH LENGTH C C COMPONENTS OF THE INDEX OF REFRACTION VECTOR C NR = N*COSEPS*COSDEL C NT = N*COSEPS*SINDEL C NP = N*SINEPS C DERY(1) = (DR/DS) = (DR/DT)*(DT/DS) = (DR/DT)*(1/VG) DERY(1)=(COSEPS*COSDEL-0.5*DN2DPS*(BO*COP*COSEPS*COSDEL-BR) 1 /(M2*BO*SIP))*RAY DR1=DERY(1) C DERY(2) = DTHETA/DS DERY(2)=(COSEPS*SINDEL-.5*DN2DPS*(BO*COP*COSEPS*SINDEL-BT) 1 /(M2*BO*SIP))*RAY/Y1 DR2=DERY(2) C DERY(3)=DPHI/DS DERY(3)=(SINEPS-0.5*DN2DPS*(BO*COP*SINEPS-BP) 1 /(M2*BO*SIP))/(Y1*SINT)*RAY DR3=DERY(3) C DERF = DNR/DS DERF=(N*RAY) * (0.5*DN2DPS*DPSIDR/M2 + 1 N*COSEPS*SINDEL*(DERY(2)/(N*RAY)) + 2 N*SINEPS*(DERY(3)/(N*RAY))*SINT) C DERG = DNTHETA/DS DERG=(N*RAY/Y1) * (0.5*DN2DPS*DPSIDT/M2 - 1 N*COSEPS*SINDEL*(DERY(1)/(N*RAY)) + 2 Y1*N*SINEPS*(DERY(3)/(N*RAY))*COST) C DERH = DNPHI/DS DERH=(N*RAY/(Y1*SINT)) * (0.5*DN2DPS*DPSIDP/M2 - 1 N*SINEPS*((DERY(1)/(N*RAY))*SINT + Y1*(DERY(2)/(N*RAY))*COST)) C DO 600 J=1,IONS1 DERF=DERF+(.5*DN2DFB(J)*DFBDR(J)+.5*DN2DFP(J)*DFP2DR(J))*RAY/N DERG=DERG+(.5*DN2DFB(J)*DFBDT(J)+.5*DN2DFP(J)*DFP2DT(J))*RAY/(Y1*N) DERH=DERH+(.5*DN2DFB(J)*DFBDP(J)+.5*DN2DFP(J)*DFP2DP(J))*RAY/ 1 (Y1*N*SINT) 600 CONTINUE C DERY(4) = DDEL/DS DERY(4)=-(SINDEL*DERF-COSDEL*DERG)/(N*COSEPS) C DERY(5) = DEPS/DS DERY(5)=-(SINEPS*(COSDEL*DERF+SINDEL*DERG)-COSEPS*DERH)/N RETURN END C C INSERT NEW MAGNET SUBROUTINE FOR JUPITER 04/28/94 C SUBROUTINE MAGNETDRIVER(DR,DTHETA,DPHI) C THIS ROUTINE CALLS THE O6 AND O6_SHEET ROUTINES AND CALCULATES THE C MAGNETIC FIELD COMPONENTS OF JUPITER INCLUDING THE CURRENT SHEET C CORRECTION (SEE CONNERNEY ET AL., , <86>, 3623, 1982. C DR IS RADIAL DISTANCE IN RJ C DTHETA AND DPHI ARE THETA AND PHI IN RADIANS C C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /BLK001/ DEGRAD,RADEG,BR,BT,BP,BO,DBODR,DBODT,DBODP, # DBRDR,DBRDT,DBRDP,DBTDR,DBTDT,DBTDP,STEPMX C C C DEFINE CONSTANTS FOR SHEET CURRENT (TAKEN FROM CONNERNEY ET AL. , C <87>, 3625, 1982. C DR0=5.0 !INNER EDGE OF CURRENT SHEET (RJ) DR1=50.0 !OUTER EDGE OF CURRENT SHEET (RJ) DD=2.5 !DISC HALF-THICKNESS (RJ) DC=.00225 !SCALE CONSTANT FOR THE CURRENT DENSITY (GAUSS) DXT=9.6 !TILT IN THETA DIRECTION OF CURRENT SHEET (DEGREES) DXP=202.0 !TILT IN PHI DIRECTION OF CURRENT SHEET (DEGREES) C CALL SHTFLD(DR,DTHETA,DPHI,DR0,DR1,DD,DC,DXT,DXP,DBR,DBT,DBP) CALL MAGNETO6(DR,DTHETA,DPHI) !CONNERNEY FIELD MODEL C C NOW ADD INTERNAL AND EXTERNAL FIELDS C BRT=BR+DBR BTT=BT+DBT BPT=BP+DBP BR=BRT BT=BTT BP=BPT BO=SQRT(BR*BR+BT*BT+BP*BP) C RETURN END C SUBROUTINE SHTFLD(DR,DTHETA,DPHI,DR0,DR1,DD,DC,DXT,DXP,DBR,DBT, * DBP) C C C C 05/06/80 REV 00 CONNERNEY C SUBROUTINE ACCEPTS DOUBLE PRECISION INPUT C AND RETURNS FIELD DUE TO SHEET FOR L-SHELL C PROGRAM COMPUTATIONS. C 05/10/80 REV 01 CHANGE TO SINGLE PRECISION FOR ADDITION TO C GENERALIZED INVERSE PROGRAM .............. C C 09/20/94 REV 02 CHANGE TO IMPLICIT REAL*8 FOR COMPATIBILITY C WITH RAY TRACING CODE C C DR : REAL*4 RADIAL POSITION IN UNITS RJ , SYS III C DTHETA: REAL*4 THETA COORDINATE IN RADIANS , USUAL SPHERICAL SYS C DPHI : REAL*4 PHI COORDINATE IN RADIANS , EAST LONGITUDE C DR0 : REAL*4 R0 , MODEL PARAMETER INNER EDGE OF SHEET C DR1 : REAL*4 R1 , MODEL PARAMETER OUTER EDGE OF SHEET C DD : REAL*4 D , MODEL PARAMETER HALF-THICKNESS OF SHEET C DC : REAL*4 CURRNT, MODEL CURRENT DENSITY C DXT : REAL*4 XT , THETA COORDINATE OF MAGNETIC DIPOLE AXIS C IN DEGREES. ALWAYS USE 9.6 DEGREES. C DXP : REAL*4 XP , PHI COORDINATE OF MAGNETIC DIPOLE AXIS C IN DEGREES SYSTEM III. 1957: 232 . 1965 : 202 ..... C DBR : REAL*4 SYSIII R,THETA,PHI COMPONENTS OF SHEET FIELD C DBT C DBP C C COMPUTES THE FIELD DUE TO THE EXTERNAL CURRENT SYSTEM C AND RETURNS THE RESULT IN GAUSS. C UTILIZES THE ANALYTICAL APPROXIMATIONS IN ALL REGIONS C AND REQUIRES SUBROUTINES : ASHORT,ALONG,VMULFF (IMSL5) C GOCART,GOCYL.......................................... C C C CM REAL DR,DTHETA,DPHI,DR0,DR1,DD,DXT,DXP,DBR,DBT,DBP CM REAL AX(3,3),AXT(3,3),B(3),S(3,3),B2(3),PHI,PV(3),PV2(3) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION AX(3,3),AXT(3,3),B(3),S(3,3),B2(3),PV(3),PV2(3) DATA TWOPIE/6.28319/ C C CONVERT TO SINGLE PRECISION C R = DR THETA = DTHETA PHI = DPHI R0 = DR0 R1 = DR1 D = DD CURRNT = DC XT = DXT XP = DXP C XT=XT*TWOPIE/360. XP=XP*TWOPIE/360. XP=TWOPIE-XP CT=COS(XT) ST=SIN(XT) CP=COS(XP) SP=SIN(XP) C C CONVERT ABOVE THE MAG. AXIS LONGITUDE TO EAST LONG. C C C BELOW WE HAVE THE TRANSFORMATION MATRIX FROM CARTESIAN C SYSTEM THREE TO CARTESIAN MAGNETIC EQUATORIAL COORDINATES. C C AX(1,1)=CT*CP AX(1,2)=SP*CT AX(1,3)=-ST AX(3,3)=CT AX(2,3)=0.0 AX(2,1)=-SP AX(2,2)=CP AX(3,1)= ST*CP AX(3,2)= ST*SP C C C BELOW WE COMPUTE THE INVERSE TRANSFORMATION ( A TRANSPOSE ) C DO 10 I=1,3 DO 11 J=1,3 AXT(I,J)=AX(J,I) 11 CONTINUE 10 CONTINUE C C C C CALL GOCART(R,THETA,PHI,PV(1),PV(2),PV(3)) CALL MATMUL(AX,PV,PV2) CALL GOCYL(PV2(1),PV2(2),PV2(3),RHO,ANGL,Z) C C ABOVE IS THE COORDINATES OF THE POINT IN MAGNETIC C EQUATORIAL COORDINATES, RHO, ANGL, Z SYSTEM. C IF(RHO .LT. R0) CALL ASHORT(RHO,Z,R0,D,AR,AZ) IF(RHO .GT. R0) CALL ALONG(RHO,Z,R0,D,AR,AZ) CALL ASHORT(RHO,Z,R1,D,AR2,AZ2) SUMR = AR-AR2 SUMZ = AZ-AZ2 C HAVE BELOW THE COMPONENTS OF SHEET B IN MAGNETIC EQUATORIAL C COORDINATE SYSTEM RHO PHI AND Z. C B(1)=COS(ANGL)*SUMR B(2)=SIN(ANGL)*SUMR B(3)=SUMZ C C CONVERT CARTESIAN MAGNETIC EQUATORIAL TO SYS 3 CARTESIAN. C CALL MATMUL(AXT,B,B2) C C ASSEMBLE THE S MATRIX FOR TRANSFORMATION FROM CARTESIAN C TO SPHERICAL COORDINATES... C CTH=COS(THETA) STH=SIN(THETA) CPH=COS(PHI) SPH=SIN(PHI) C S(1,1) = STH*CPH S(1,2) = STH*SPH S(1,3) = CTH S(2,1) = CTH*CPH S(2,2) = CTH*SPH S(2,3) = - STH S(3,1) = - SPH S(3,2) = CPH S(3,3) = 0.0 C C CONVERT SYS 3 CARTESIAN TO SYS 3 SPHERICAL COORDINATES. C CALL MATMUL(S,B2,B) C C C NOW SCALE ALL B UP TO GAMMAS WITH UI SCALE FACTOR. C SBR=B(1)* CURRNT SBT=B(2)* CURRNT SBP=B(3)* CURRNT C DBR = SBR * 1. DBT = SBT * 1. DBP = SBP * 1. C RETURN END SUBROUTINE GOCART(R,THETA,PHI,X,Y,Z) C Converts spherical position vector to Cartesian. C This subroutine was included with Dr. Connerney's SHTFLD package. CM REAL PHI C IMPLICIT REAL*8 (A-H,O-Z) C X=R*SIN(THETA)*COS(PHI) Y=R*SIN(THETA)*SIN(PHI) Z=R*COS(THETA) RETURN END SUBROUTINE MATMUL(A,B,C) C Multiplies a 3x3 matrix A by a 3-vector B. CM REAL A(3,3),B(3),C(3) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION A(3,3),B(3),C(3) C(1)=0. C(2)=0. C(3)=0. DO 10 I=1,3 DO 20 K=1,3 C(I)=C(I)+A(I,K)*B(K) 20 CONTINUE 10 CONTINUE RETURN END SUBROUTINE GOCYL(X,Y,Z,RHO,PHI,ZED) C Converts Cartesian position vector to cylindrical. C This subroutine was included with Dr. Connerney's SHTFLD package. CM REAL PHI C IMPLICIT REAL*8 (A-H,O-Z) C DATA TWOPIE/6.28319/ RHO=SQRT( X**2 + Y**2 ) PHI=ATAN2(Y,X) IF (PHI.LT.0) PHI=TWOPIE + PHI ZED=Z RETURN END C SHORT : PROG. TO COMPUTE SHEET FIELDS NEAR AXIS C C 04/04/80 REV 00 C This subroutine was included in Dr. Connerney's SHTFLD package. C SUBROUTINE ASHORT(RHO,Z,A,D,BR,BZ) C IMPLICIT REAL*8 (A-H,O-Z) C C F1 = ((Z-D)**2 + A**2 ) F1 = SQRT(F1) F2 = ((Z+D)**2 + A**2 ) F2 = SQRT(F2) F3 = SQRT( Z**2 + A**2 ) C BR = RHO*(1/F1 - 1/F2 ) / 2 BZ = 2*D/F3 - RHO**2 * ( (Z-D)/F1**3-(Z+D)/F2**3)/4 C RETURN END C LONG : PROG TO COMPUTE SHEET FIELDS IN THE REGION C RHO>>R0, I.E., THE FAR AWAY FIELDS............... C C This subroutine was included in Dr. Connerney's SHTFLD package. C 4/5/80 REV 00 C 07/07/81 REV 01 : CORRECT B RHO FOR Z MAG > D, Z NEG C SUBROUTINE ALONG(RHO,Z,A,D,BR,BZ) C IMPLICIT REAL*8 (A-H,O-Z) C C F1 = ((Z-D)**2 + RHO**2 ) F2 = ((Z+D)**2 + RHO**2 ) C F1 = SQRT( F1 ) F2 = SQRT( F2 ) C BR = ( F1 - F2 + 2*D )/ RHO IF ( ABS(Z) .LT. D ) BR = ( F1 - F2 + 2*Z ) / RHO IF ( ABS(Z) .GT. D .AND. Z .LT. 0. ) BR=(F1-F2-2*D)/RHO C BR = BR -( A**2 * RHO / 4)*(1/F1**3 - 1/F2**3) C BZ = 2 * D / SQRT( Z**2 + RHO**2 ) BZ = BZ - ( A**2/4 )*( (Z-D)/F1**3 - (Z+D)/F2**3 ) C RETURN END C SUBROUTINE O4RFT(R,T,P,BR,BT,BP) C C C SPHERICAL HARMONIC FUNCTION EXPANDED ABOUT (X0,Y0,Z0) C AND GRADIENT EVALUATED AT (R,T,P) C C R = RADIAL DISTANCE ORIGIN=(0,0,0) C T = COLATITUDE (THETA) (RADIANS) ORIGIN=(0,0,0) C P = EAST LONGITUDE (PHI) (RADIANS) ORIGIN=(0,0,0) C C GRADIENT COMPONENTS BR,BT,BP ARE SPHERICAL POLAR COMPONENTS C WITH RESPECT TO ORIGIN AT (X0,Y0,Z0) C C INT = DEGREE OF INTERNAL SOURCE, COEFFICIENTS GIJ, HIJ C EXT = DEGREE OF EXTERNAL SOURCE, COEFFICIENTS EIJ, FIJ C C MODIFICATIONS 13 JAN 75 C 1 AVOID TRANSLATION CODE WHEN X0=Y0=Z0=APPROX ZERO C 2 AVOID Y/RO, X/RO WHEN X=Y=RO=APPROX ZERO C C MODIFICATIONS 5 MAR 75 C 1 CODE FOR 4TH DEGREE INTERNAL C 2 REORDER COMPONENT ACCUMULATION, SMALLEST TO LARGEST (R.GT.1) C I4,I3,I2,I1,E1,E2,..... C C MODIFICATIONS 23 OCT 85 C 1 REMOVE COMMON BLOCK, ASSIGN PARAMETERS INTERNALLY C C C IMPLICIT REAL*8 (A-H,O-Z) C INTEGER EXT DATA SQ3/1.73205/,SQ5/2.23607/,SQ8/2.82843/,SQ15/3.87298/ DATA SQ7/2.64575/ INT = 3 EXT = 0 X0=0.0 Y0=0.0 Z0=0.0 G10 = 4.218 G11 =-0.664 H11 = 0.264 G20 =-0.203 G21 =-0.735 H21 =-0.469 G22 = 0.513 H22 = 0.088 G30 =-0.233 G31 =-0.076 H31 =-0.580 G32 = 0.168 H32 = 0.487 G33 =-0.231 H33 =-0.294 C G40 = 0.0 G41 = 0.0 H41 = 0.0 G42 = 0.0 H42 = 0.0 G43 = 0.0 H43 = 0.0 G44 = 0.0 H44 = 0.0 E10 = 0.0 E11 = 0.0 F11 = 0.0 E20 = 0.0 E21 = 0.0 F21 = 0.0 E22 = 0.0 F22 = 0.0 C C ST = SIN(T) CT = COS(T) SP = SIN(P) CP = COS(P) RR = R EPS = 1.0E-15 IF(ABS(X0) - EPS) 1,1,3 1 IF(ABS(Y0) - EPS) 2,2,3 2 IF(ABS(Z0) - EPS) 18,18,3 C C COMPUTE CARTESIAN COORDINATES OF (R,T,P) IN SYSTEM HAVING C ORIGIN AT (X0,Y0,Z0) 3 CONTINUE X = R*ST*CP - X0 Y = R*ST*SP - Y0 Z = R*CT - Z0 IF(ABS(X) - EPS) 4,4,6 4 IF(ABS(Y) - EPS) 5,5,6 5 ST = 0.0 CT = 1.0 SP = 0.0 CP = 1.0 GO TO 18 C C COMPUTE SPHERICAL COORDINATES OF (R,T,P) C IN SYSTEM HAVING ORIGIN AT (X0,Y0,Z0) 6 CONTINUE RO = SQRT(X*X + Y*Y) RR = SQRT(RO*RO + Z*Z) ST = RO/RR SIN(THET CT = Z/RR COS(THET SP = Y/RO SIN(PHI) CP = X/RO COS(PHI) 18 CONTINUE C COMPUTE ASSOCIATED LEGENDRE POLYNOMIALS ETC. C SCHMIDT QUASI-NORMALIZED, DEGREE 1, 2, 3, 4 AS NEEDED. C C 1ST AND 2ND DEGREE FUNCTIONS R3 = 1.0/RR**3 R4 = R3/RR S2P = 2.0*SP*CP C2P = 2.0*CP*CP - 1.0 ST2 = ST*ST CT2 = CT*CT P20 = 1.5*CT2 - 0.5 P21 = SQ3*ST*CT P22 = 0.5*SQ3*ST2 DP20 =-3.0*ST*CT DP21 = SQ3*(CT2 - ST2) DP22 = SQ3*ST*CT C REMOVE SIN FACTOR FOR PHI COMPONENTS P21S = SQ3*CT P22S = 0.5*SQ3*ST C IF(INT-3) 7,8,8 7 IF(EXT-3) 9,8,8 C C 3RD DEGREE FUNCTIONS C 8 R5 = R4/RR P30 = CT*(2.5*CT2 - 1.5) P31 = SQ3*ST*(5.0*CT2 - 1.0)/SQ8 P32 = 0.5*SQ15*CT*ST2 P33 = SQ5*ST*ST2/SQ8 DP30 = ST*(-7.5*CT2 + 1.5) DP31 = SQ3*CT*(5.0*CT2 - 10.0*ST2 - 1.0)/SQ8 DP32 = SQ15*ST*(CT2 - (ST2/2.0)) DP33 = 3.0*SQ5*ST2*CT/SQ8 C REMOVE SIN FACTOR FOR PHI COMPONENTS P31S = SQ3*(5.0*CT2 - 1.0)/SQ8 P32S = 0.5*SQ15*ST*CT P33S = SQ5*ST2/SQ8 S3P = SP*C2P + CP*S2P C3P = CP*C2P - SP*S2P 9 CONTINUE IF(INT - 4) 19,20,20 19 IF(EXT - 4) 21,20,20 C C 4TH DEGREE FUNCTIONS 20 CONTINUE R6 = R5/RR S4P = SP*C3P + CP*S3P SIN(4P) C4P = CP*C3P - SP*S3P COS(4P) P40 = 0.125*(3.0 + CT2*(35.0*CT2 - 30.0)) P41 = (7.0*CT*P31 - SQ8*P21)/SQ15 P42 = (7.0*CT*P32 - SQ5*P22)/(2.0*SQ3) P43 = SQ7*CT*P33 C *** P44 = SQ7 *ST*P33/SQ8 P44 = 0.935414*ST*P33 DP40 = 0.25*ST*CT*(-70.0*CT2 + 30.0) DP41 = (7.0*(CT*DP31 - ST*P31) - SQ8*DP21)/SQ15 DP42 = (7.0*(CT*DP32 - ST*P32) - SQ5*DP22)/(2.0*SQ3) DP43 = SQ7*(CT*DP33 - ST*P33) C *** DP44 = SQ35*CT*ST*ST2 / 2.0 DP44 = 2.95804*CT*ST*ST2 C REMOVE SIN FACTOR FOR PHI COMPONENTS P41S = (7.0*CT*P31S - SQ8*P21S)/SQ15 P42S = (7.0*CT*P32S - SQ5*P22S)/(2.0*SQ3) P43S = SQ7*CT*P33S P44S = 0.935414*ST*P33S 21 CONTINUE C EVALUATE GRADIENT OF SPHERICAL FUNCTION AT (R,T,P) C SPHERICAL COMPONENTS BR,BT,BP ARE WITH RESPECT TO C COORDINATE ORIGIN AT (X0,Y0,Z0) BR = 0.0 BT = 0.0 BP = 0.0 IF(INT) 13,13,22 22 IF(INT - 4) 24,23,23 C INTERNAL 4TH DEGREE 23 BR=BR+ 5.0*R6*(G40*P40 + (G41* CP + H41* SP)*P41 C + (G42*C2P + H42*S2P)*P42 C + (G43*C3P + H43*S3P)*P43 C + (G44*C4P + H44*S4P)*P44 ) C BT=BT- R6*(G40*DP40+ (G41* CP + H41* SP)*DP41 C + (G42*C2P + H42*S2P)*DP42 C + (G43*C3P + H43*S3P)*DP43 C + (G44*C4P + H44*S4P)*DP44) C BP=BP+ R6*( (G41* SP - H41* CP)*P41S C +2.0*(G42*S2P - H42*C2P)*P42S C +3.0*(G43*S3P - H43*C3P)*P43S C +4.0*(G44*S4P - H44*C4P)*P44S) C 24 IF(INT - 3) 25,12,12 C INTERNAL 3RD DEGREE 12 BR=BR+ 4.0*R5*(G30*P30 + (G31* CP + H31* SP)*P31 C + (G32*C2P + H32*S2P)*P32 C + (G33*C3P + H33*S3P)*P33 ) C BT=BT- R5*(G30*DP30+ (G31* CP + H31* SP)*DP31 C + (G32*C2P + H32*S2P)*DP32 C + (G33*C3P + H33*S3P)*DP33 ) C BP=BP+ R5*( (G31* SP - H31* CP)*P31S C +2.0*(G32*S2P - H32*C2P)*P32S C +3.0*(G33*S3P - H33*C3P)*P33S ) C 25 IF(INT - 2) 10,11,11 C INTERNAL 2ND DEGREE 11 BR=BR +3.0*R4*(G20*P20 + (G21* CP + H21* SP)*P21 C + (G22*C2P + H22*S2P)*P22) C BT=BT - R4*(G20*DP20+ (G21* CP + H21* SP)*DP21 C + (G22*C2P + H22*S2P)*DP22) C BP=BP + R4*( (G21* SP - H21* CP)*P21S C +2.0*(G22*S2P - H22*C2P)*P22S) C C INTERNAL 1ST DEGREE 10 BR=BR+ 2.0*R3*(G10* CT + (G11* CP + H11* SP)*ST ) BT=BT- R3*(-G10*ST + (G11* CP + H11* SP)*CT ) BP=BP+ R3*( G11* SP - H11* CP ) 13 IF(EXT) 17,17,14 C EXTERNAL 1ST DEGREE 14 BR=BR -E10* CT - (E11* CP + F11* SP)*ST BT=BT +E10* ST - (E11* CP + F11* SP)*CT BP=BP + E11* SP - F11* CP IF(EXT-1) 17,17,15 C EXTERNAL 2ND DEGREE 15 BR=BR -2.0*RR*(E20*P20 + (E21* CP + F21* SP)*P21 C + (E22*C2P + F22*S2P)*P22) C BT=BT -RR*(E20*DP20+ (E21* CP + F21* SP)*DP21 C + (E22*C2P + F22*S2P)*DP22) C BP=BP + RR*( (E21* SP - F21* CP)*P21S C +2.0*(E22*S2P - F22*C2P)*P22S) IF(EXT-2) 17,17,16 C EXTERNAL 3RD DEGREE 16 CONTINUE C INSERT CODE FOR EXT = 3 17 CONTINUE RETURN END C################################################################### C INSERT NEW MAGNET SUBROUTINE FOR JUPITER 04/27/94 C C Y1 IS RADIAL DISTANCE IN RJ C Y2, Y3 ARE THETA AND PHI IN RADIANS C SUBROUTINE MAGNETO6(Y1,Y2,Y3) !CONNERNEY FIELD MODEL IMPLICIT REAL*8 (A-H,O-Z) !### THIS LINE NEEDED FOR OTHER CODES #### CM DIMENSION Y(6) COMMON /BLK001/ DEGRAD,RADEG,BR,BT,BP,BO,DBODR,DBODT,DBODP, # DBRDR,DBRDT,DBRDP,DBTDR,DBTDT,DBTDP,STEPMX COMMON/BLKDIP/IDIP,THETADIP,PHIDIP,THETMAG,BLD,IMAG C COMMON/LSHDEN/EL0,CIL,EKLP C COMMON/OTDBLK/LSHELL,DLDR,DLDT DIMENSION GG(9),HH(9),PP(9),DP(9),AOR(4),SP(4),CP(4),DP2(9) C C O-6 MODEL COEFF. C GG(1)=4.24202 GG(2)=-0.65929 GG(3)=-0.02181 GG(4)=-0.71106 GG(5)=0.48714 GG(6)=0.07565 GG(7)=-0.15493 GG(8)=0.19775 GG(9)=-0.17958 HH(1)=0.0 HH(2)=0.24116 HH(3)=0.0 HH(4)=-0.40304 HH(5)=0.07179 HH(6)=0.0 HH(7)=-0.38824 HH(8)=0.34243 HH(9)=-0.22439 C PRAD=71300. !RADIUS IN KM C C MAGNETIC COEFFICIENTS FOR THE OLD O4 MODEL CM IF(IMAG.EQ.2) THEN GG(1)=4.22 GG(2)=-.442 GG(3)=-.203 GG(4)=-.871 GG(5)=.331 GG(6)=-.233 GG(7)=-.357 GG(8)=.506 GG(9)=-.292 HH(1)=0.0 HH(2)=0.562 HH(3)=0.0 HH(4)=-.037 HH(5)=-.402 HH(6)=0.0 HH(7)=-.463 HH(8)=.096 HH(9)=.233 ENDIF !IMAG=2 C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C MAGNETIC MODEL COEFF. FOR A DIPOLE OR O6 C C DEFINE ROTATED AND TILTED COORDINATES WRT SYSTEM III COORDINATES C C SYSIII COORDS: C X3D=Y1*SIN(Y2)*COS(Y3) Y3D=Y1*SIN(Y2)*SIN(Y3) Z3D=Y1*COS(Y2) C THETADIP=9.6 !DIPOLE DIP PHIDIP=160.0 !DIPOLE AZIMUTH ANGLE C CTD=COSD(THETADIP) STD=SIND(THETADIP) CPD=COSD(PHIDIP) SPD=SIND(PHIDIP) C C TILTED AND ROTATED COORDS. C CM XDIP=CPD*X3D+SPD*Y3D CM YDIP=-SPD*CTD*X3D+CTD*CPD*Y3D+STD*Z3D CM ZDIP=STD*SPD*X3D-STD*CPD*Y3D+CTD*Z3D C XDIP=CTD*CPD*X3D+CTD*SPD*Y3D-STD*Z3D YDIP=-SPD*X3D+CPD*Y3D ZDIP=STD*CPD*X3D+STD*SPD*Y3D+CTD*Z3D C DIPFAC=XDIP*XDIP+YDIP*YDIP ROWDIP=SQRT(DIPFAC) RDIP=SQRT(DIPFAC+ZDIP*ZDIP) STDIP=ROWDIP/RDIP THETMAG=ASIND(STDIP) BLD=Y1/(STDIP*STDIP) !L-VALUE ONLY FOR MAGNETIC COORDS. C C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C RS=Y1 ! RS IS THE DISTANCE IN RJ Ss=Y2 Ss=SIN(Ss) IF(ABS(Ss).LT.1.0D-12) THEN IF(Ss.LT.0.0) THEN Ss=-1.0D-12 ELSE Ss= 1.0D-12 END IF END IF C=Y2 C=COS(C) PHI=Y3 ! PHI IS DEFINED IN A RIGHT HANDED COORDINATE SYSTEM DO 5 I=1,9 ! ZERO OUT ARRAYS DP(I)=0.0 PP(I)=0.0 5 CONTINUE DO 6 I=1,4 SP(I)=0.0 CP(I)=0.0 6 CONTINUE DO 7 I=1,3 AOR(I)=0.0 7 CONTINUE x=2.0 X=SQRT(x) z=3.0 Z=SQRT(z) ys=5.0 Ys=SQRT(ys) IF(RS.EQ.0.0) RS=1.0E-03 AR=1.0/RS SP(1)=0.0 CP(1)=1.0 SP(2)=SIN(PHI) CP(2)=COS(PHI) AOR(1)=AR**3 ! AOR IS THE RADIAL DISTANCE ARRAY CUBED,FORTH,FIFTH AOR(2)=AOR(1)*AR AOR(3)=AOR(2)*AR AOR(4)=AOR(3)*AR DO 90 M=3,4 N=M-1 SP(M)=SP(2)*CP(N)+CP(2)*SP(N) CP(M)=CP(2)*CP(N)-SP(2)*SP(N) 90 CONTINUE B1=0.0 B2=0.0 B3=0.0 BTDR=0.0 BTDT=0.0 BTDP=0.0 BRDR=0.0 BRDT=0.0 BRDP=0.0 BPDR=0.0 BPDT=0.0 BPDP=0.0 PP(1)=C ! LEGENDRE POLY AND THEIR DERIVATIVES PP(2)=Ss PP(3)=(3.0*C*C-1.0)/2.0 PP(4)=Z*C*Ss PP(5)=Z*Ss*Ss/2.0 PP(6)=C*(5.0*C*C-3.0)/2.0 PP(7)=Z*Ss*(5.0*C*C-1.0)/(2.0*X) PP(8)=Ys*Z*C*Ss*Ss/2.0 PP(9)=Ys*Ss*Ss*Ss/(2.0*X) DP(1)=1.0 DP(2)=-C/Ss DP(3)=3.0*C DP(4)=Z*(1.0-2.0*C*C)/Ss DP(5)=-Z*C DP(6)=(15.0*C*C-3.0)/2.0 DP(7)=Z*C*(11.0-15.0*C*C)/(2.0*X*Ss) DP(8)=Ys*Z*(1.0-3.0*C*C)/2.0 DP(9)=-3.0*Ys*C*Ss/(2.0*X) DP2(1)=0.0 DP2(2)=-1.0/(Ss*Ss*Ss) DP2(3)=3.0 DP2(4)=Z*C*(2.0*C*C-3.0)/(Ss*Ss*Ss) DP2(5)=-Z DP2(6)=15.0*C DP2(7)=Z*(11.0-45.0*C*C+30.0*C*C*C*C)/(2.0*X*Ss*Ss*Ss) DP2(8)=-3.0*Ys*Z*C DP2(9)=-3.0*Ys*(1.0-2.0*C*C)/(2.0*X*Ss) DO 32 N=1,3 FN=N+1 FN2=N+2 SUMR=0.0 SUMT=0.0 SUMP=0.0 SUMT3=0.0 SUMT4=0.0 SUMP4=0.0 DO 33 M=1,4 K=M-1 FM=K IF (K.GT.N) GO TO 33 I=(N-1)*3+K IF (I.EQ. 1) I=2 IF (I.EQ. 0) I=1 TS=GG(I)*CP(M)+HH(I)*SP(M) TS2=-GG(I)*SP(M)+HH(I)*CP(M) SUMR=SUMR+PP(I)*TS SUMT=SUMT+DP(I)*TS SUMP=SUMP+FM*PP(I)*TS2 SUMT3=SUMT3+C*(DP(I)*TS)-Ss*Ss*DP2(I)*TS SUMT4=SUMT4+DP(I)*TS2*FM SUMP4=SUMP4+FM*FM*PP(I)*TS 33 CONTINUE B1=B1+AOR(N)*FN*SUMR ! B FIELD DUE TO INTERMAL SOURCES ONLY B2=B2+AOR(N)*SUMT B3=B3-AOR(N)*SUMP BTDR=BTDR+FN2*AOR(N+1)*SUMT BTDT=BTDT+AOR(N)*SUMT3 BTDP=BTDP+AOR(N)*SUMT4 BRDR=BRDR+FN*FN2*AOR(N+1)*SUMR BRDT=BRDT+FN*AOR(N)*SUMT BRDP=BRDP+FN*AOR(N)*SUMP BPDR=BPDR+FN2*AOR(N+1)*SUMP BPDT=BPDT+AOR(N)*SUMT4 BPDP=BPDP+AOR(N)*SUMP4 32 CONTINUE BP=B3/Ss BT=B2*Ss BR=B1 BTDR=-BTDR*Ss/PRAD BTDP=BTDP*Ss DBRDR=-BRDR/PRAD DBRDT=-Ss*BRDT DBRDP=BRDP BPDR=BPDR/(Ss*PRAD) BPDT=BP*C/Ss+BPDT BPDP=BPDP/Ss BO=SQRT(BR*BR+BT*BT+BP*BP) DBODR=(BR*DBRDR+BT*BTDR+BP*BPDR)/BO DBODT=(BR*DBRDT+BT*BTDT+BP*BPDT)/BO DBODP=(BR*DBRDP+BT*BTDP+BP*BPDP)/BO DBTDR=BTDR DBTDT=BTDT DBTDP=BTDP RETURN END C C @@@@@@@@@@@@@@ INCLUDE NEW DIVINE AND GARRETT PLASMA MODEL @@@@ C CM SUBROUTINE MODEL(Y,FP2,DFP2DR,DFP2DT,DFP2DP,FB,FB2,DFBDR,DFBDT, CM #DFBDP,ANER) SUBROUTINE MODEL(Y) INCLUDE 'JUPCOMMONS.FOR' C CM IMPLICIT REAL*8 (A-H,O-Z) C CM REAL*8 NR,NT,NP CM DIMENSION Y(6),FP2(8),DFP2DR(8),DFP2DT(8),FB(8),FB2(8), CM #DFBDR(8),DFBDT(8),PLASCON(8),ZI(8) CM DIMENSION DFP2DP(8),DFBDP(8),W(8) CM DIMENSION DELFBR(8),DELFBT(8),DELFBP(8) CM COMMON /INDEX/CN,S,P,D,SIP,COP,ALTUP,ALTDN,F,ALT,ACC, CM 1 COSDEL,SINDEL,RO,W,NR,NT,NP,FP,FG,T,ALT1, CM 2 ALT2,IC,IONS,IPRINT,ITURN,ITMAX,ILINES,IT,PRAD,JK,ISERVO,STEP CM COMMON /BLK001/ DEGRAD,RADEG,BR,BT,BP,BO,DBODR,DBODT,DBODP, CM # DBRDR,DBRDT,DBRDP,DBTDR,DBTDT,DBTDP,STEPMX CM COMMON/BLKA1/ZI !MAIN AND MODEL CM COMMON/DGPLAS/PLASCON,IONS1 CCCM COMMON/BLKOMNI/IFLAGG,IRADEX !IN ALL ROUTINES C C DEFINE THE PLASMA CONSTANT [FP2=PLASCON*Ne] C IONS1=IONS+1 C C DEFINE RADIAL DISTANCE IN RJ C RJ=Y(1)/PRAD C PI=3.1415927 DO JS=1,IONS1 IF(JS.EQ.1) THEN PLASCON(JS)=8.06404E+07 ELSE !JS NE 1 PLASCON(JS)=4.41E+04*ZI(JS-1)*ZI(JS-1)/W(JS-1) !ZI IS CHARGE STATE ENDIF !JS ENDDO !JS C CONS(1)=8.0503E 07 C CONS(2)=4.3915E 04/(W(1)) C CONS(3)=4.3915E 04/(W(2)) C CONS(4)=4.3915E 04/W(3) C C WRITE(6,*) ' RJ,Y2,Y3 ',RJ,Y(2),Y(3) CM CALL MAGNET(RJ,Y(2),Y(3)) CALL MAGNETDRIVER(RJ,Y(2),Y(3)) C WRITE(6,*) ' BO= ',BO C C ELECTRON AND ION GYROFREQUENCIES C C WRITE(6,3332) IONS,IONS1 C 3332 FORMAT(1X,'IONS,IONS1',2I3) DO J=1,IONS1 IF(J.EQ.1) THEN CON1=2.7952E06 CON2=1.0 ELSE CON1=1.5224E03 CON2=W(J-1) ENDIF !J=1 C WRITE(6,3333) J,CON1,CON2 C 3333 FORMAT(1X,'J,CON1,CON2',I4,1P2E12.3) FB(J)=ZI(J)*BO*CON1/CON2 !ZI IS CHARGE STATE FB2(J)=FB(J)*FB(J) DFBDR(J)=ZI(J)*DBODR*CON1/CON2 DFBDT(J)=ZI(J)*DBODT*CON1/CON2 DFBDP(J)=ZI(J)*DBODP*CON1/CON2 ENDDO !J LOOP WLONG=360.0-Y(3)*DEGRAD ALAT=(PI/2.0-Y(2))*DEGRAD CALL LOPLAS(RJ,WLONG,ALAT,TEMP,RHOT,RHOE,RHOH1, #RH001,RH002,RHOS1,RHOS2,RHOS3,RHONA,VCNVC,DFP2DR,DFP2DT, #DFP2DP,FP2) FP=SQRT(FP2(1)) C C INSERT TO TEST DERIVATIVES C CC DELY1=ABS(Y(1)-OLDY1) CC DELY2=ABS(Y(2)-OLDY2) CC DELY3=ABS(Y(3)-OLDY3) CC DELFP=ABS(OLDFP2-FP2(1)) CC IF(DELY1.NE.0.0) TESTDER1=DELFP/DELY1 CC IF(DELY2.NE.0.0) TESTDER2=DELFP/DELY2 CC IF(DELY3.NE.0.0) TESTDER3=DELFP/DELY3 CC WRITE(6,1301) TESTDER1,TESTDER2,TESTDER3 CC 1301 FORMAT(1X,'TEST DFP2/DR,/DT/DP= ',1P3E12.4) CC WRITE(6,1303) DFP2DR(1),DFP2DT(1),DFP2DP(1) CC 1303 FORMAT(1X,'DFP2DR/DT/DP= ',1P3E12.4) CC OLDFP2=FP2(1) CC OLDY1=Y(1) CC OLDY2=Y(2) CC OLDY3=Y(3) C C CHANGE TO "OUTP" NOTATION C ANER=RHOE RETURN END C C HERE ARE THE SUBROUTINES C SUBROUTINE LOPLAS(RJ,WLONG,ALAT,TEMP,RHOT,RHOE,RHOH1, #RH001,RH002,RHOS1,RHOS2,RHOS3,RHONA,VCNVC, #DFP2DR,DFP2DT,DFP2DP,FP2) CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8) C C WRITTEN BY NEIL DEVINE AND HANK GARRETT FROM JGR,88,6889,1983 C C RJ = Radial distance position in Rj C WLONG = System III Jovian graphic longitude C ALAT = System III Jovian graphic latitude C RETURNED ARGUMENTS C TEMP = Temperature of the plasma C RHOT = Total density at the EQUATOR! C RHOE = Electron density at your position C RHOH1 = Proton density at your position C RHOO1 = Oxygen O+ density at your position C RHOO2 = O++ density C RHOS1 = Sulfer S+ density C RHOS2 = S++ density C RHOS3 = S+++ density C RHONA = Sodium Na+ density C VCNVC = Convection velocity at your position C DIMENSION R(8) C RAD=.01745 VLAT=ALAT*RAD CALL DENT(RJ,RHOT,TEMP) IF(RJ.LT.3.8)CALL INPLS(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) IF(RJ.LT.3.8)GO TO 99 IF(RJ.LT.5.5)CALL CLTOR(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) IF(RJ.LT.5.5)GO TO 99 IF(RJ.LT.7.9)CALL WMTOR(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) IF(RJ.LT.7.9)GO TO 99 IF(RJ.LT.20.)CALL INDSC(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) IF(RJ.LT.20.)GO TO 99 CALL OTDSC(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) 99 VCNVC=V RHOE=R(1) RHOH1=R(2) RH001=R(3) RH002=R(4) RHOS1=R(5) RHOS2=R(6) RHOS3=R(7) RHONA=R(8) C TYPE *,TEMP,RHOT,RHOE,RHOH1,RH001,RH002,RHOS1,RHOS2,RHOS3,RHONA RETURN END SUBROUTINE DENT (RJ,RHOT,TEMP) CM IMPLICIT REAL*16 (A-H,O-Z) CM REAL*16 LOGT(16),LOGN(16) C IMPLICIT REAL*8 (A-H,O-Z) C CM REAL LOGT(16),LOGN(16) REAL*8 LOGT(16),LOGN(16) DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8) DIMENSION RR(16) C TABLE 3 LINEAR INTERPOLATION OF LOGS DATA RR/3.8,4.9,5.1,5.3,5.5,5.65,5.8,5.9,6.4, 17.4,7.9,10.0,20.0,60.0,100.0,170.0/ DATA LOGT/1.67,-.31,-.18,.37,.92,1.15,1.33, 21.54,1.63,1.67,1.75,2.0,2.0,2.0,2.0,2.0/ DATA LOGN/1.55,2.75,2.91,3.27,2.88,3.57,3.31,3.35, 33.18,2.78,2.25,1.48,.2,-2.,-2.,-3./ RHOT=4.65 !DIVINE AND GARRETT VALUE FOR INPLS CM RHOT=3.65 !TEST VALUE ######## CM RHOT=2.325 CM RHOT=5.65 TEMP=46. IF(RJ.LT.3.8)GO TO 99 DO 1 I=1,16 IF(RJ.LT.RR(I))GO TO 2 1 CONTINUE 2 DR=(RJ-RR(I-1))/(RR(I)-RR(I-1)) ALN=LOGN(I-1)+(LOGN(I)-LOGN(I-1))*DR RHOT=10.**ALN ALT=LOGT(I-1)+(LOGT(I)-LOGT(I-1))*DR TEMP=10.**ALT 99 RETURN END SUBROUTINE LPLAS2(RJ,ALAT,WLONG,RHOEW,RHOHC,RHOHW,TEW,THC,THW, #DFP2DR,DFP2DT,DFP2DP,FP2) CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C C COMMON/DGPLAS/DFP2DR,DFP2DT,DFP2DP,FP2,PLASCON DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8) D=RJ IF (D.LT.10.) D=10. TEW =1000. CALL DENT(RJ,RHOT,TEMP) THC=TEMP THW=30000. R=10.**(EXP((-D+30.78)/16.9)-3.0) RDD=.0174532 AP=10.77*RDD IF(RJ.LT.20.)Z=RJ*TAN(AP)*COS(RDD*(WLONG-21.)) IF(RJ.GE.20.) Z=20.*TAN(AP)*COS(RDD*(WLONG-21.)-.015708*(RJ-20.)) Z0=RJ*SIN(ALAT*RDD) EDZ=0 RT=-ABS(Z-Z0)/2. IF(RT.LT.-80.) GO TO 2 EDZ=EXP(RT) 2 R=R*EDZ RHOEW=3.*R RHOHC=2.*R RHOHW=R RETURN END SUBROUTINE INPLS(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) C C THIS ROUTINE IS THE ONLY ONE NEEDING THE IONOSPHERE C DENSITY VALUES. C CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON/DGPLAS/PLASCON,IONS1 COMMON/BLKTOR/CLTORFAC,WMTORFAC,PLASFAC DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8),PLASCON(8) DIMENSION R(8),G(8) DATA G/1.0,0.0,.2,.02,.7,.03,0.0,0.0/ PRAD=71300.0 V=12.6*RJ Y1=RJ*PRAD RAD=.01745 RHI=1.D25 RLO=1.D-25 VLATC=.123*COS(RAD*(WLONG-21.)) CM EXPP=0 EXPP=RLO RT=DLOG10(RHOT)+.4342935*((7.68/RJ)-(RJ-1.)**2*(VLAT 1-VLATC)**2) IF(RT.LT.-35.) GO TO 2 IF(RT.EQ.0.0) GO TO 2 !ADDED BY JDM ON 2/27/98 EXPP=10.**RT CM 2 DO 1 I=1,8 2 DO 1 I=1,IONS1 CM 1 R(I)=G(I)*EXPP R(I)=G(I)*EXPP R(I)=PLASFAC*R(I) !ADDED BY JDM ON AUG 2, 1995 IF((R(I)).GT.(RHI)) R(I)=RHI IF((R(I)).LT.(RLO)) R(I)=RLO 1 CONTINUE C C EVALUATE THE NECESSARY DERIVATIVES C RJ1=RJ-1.0 DELLAM=VLAT-VLATC C C DEFINE IONOSPHERIC DENSITY C PIY2=3.1415927/2.0 THETION=PIBY2-VLAT HITE=(RJ-1.0)*PRAD CMM CALL IONOS(HITE,THETION,DNIONDR,FP2EION,WIONOS,WMAGOS, CMM #DWIDR,DWMDR) C CM DO JE=1,8 DO JE=1,IONS1 C DNDR=-(7.68/(RJ*RJ)+2.0*RJ1*DELLAM**2)*R(JE) DNDR=-((7.68*PRAD)/(Y1*Y1)+2.0*RJ1*DELLAM**2/PRAD)*R(JE) DNDT=2.0*DELLAM*RJ1*RJ1*R(JE) DNDPHI=2.0*RJ1*DELLAM*(.123)*SIND(WLONG-21.0)*R(JE) FP2(JE)=PLASCON(JE)*R(JE) DFP2DR(JE)=PLASCON(JE)*DNDR DFP2DT(JE)=PLASCON(JE)*DNDT DFP2DP(JE)=PLASCON(JE)*DNDPHI CM WRITE(6,54) JE,PLASCON(JE),R(JE) CM 54 FORMAT(2X,'INPLS: J,PC,R',I2,1P2E12.3) END DO !JE C C FOR IONOSPHERE WE MUST WEIGHT THE VALUES C CMM DFP2DR(1)=DWIDR*FP2EION+WIONOS*DNIONDR*PLASCON(1)+ CMM #DWMDR*FP2(1)+WMAGOS*DFP2DR(1) C C WRITE STATEMENTS C CX WRITE(6,500) DWIDR,FP2EION,WIONOS,DNIONDR,DWMDR,WMAGOS, CX # DFP2DR(1) CX 500 FORMAT(1X,'DWIDR,FP2EDR,WI,DNIDR,DWMDR,WM,DFP2DR ',/,1X, CX #1P7E12.3) C C FOLLOWING IS FOR H+ C C EL=RJ/(COS(VLAT)**2) C Z0=(7.*RJ-26.)/30.*COS((WLONG-21.)*RAD) C ENP=(30.78-10.)/16.9*(10.**-3) C R(2)=ENP*EXP(-ABS((RJ*VLAT-Z0)/2.)) C TYPE*,'L ',EL,'NP ',ENP RETURN END C SUBROUTINE CLTOR(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON/DGPLAS/PLASCON,IONS1 COMMON/BLKTOR/CLTORFAC,WMTORFAC,PLASFAC DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8),PLASCON(8) DIMENSION R(8),G(8) DATA G/1.0,0.0,.2,.02,.7,.03,0.0,0.0/ RAD=.01745 PRAD=71300.0 Y1=RJ*PRAD H=.2*SQRT(TEMP) V=12.6*RJ Z=RJ*.123*COS(RAD*(WLONG-21.)) EXPP=0 RT=-((RJ*VLAT-Z)/H)**2 IF(RT.GT.-80.) THEN EXPP=EXP(RT) ELSE WRITE(6,*) ' RT < -80 IN "CLTOR"!!' WRITE(6,33) RJ,VLAT,WLONG,Z,TEMP WRITE(6,34) G(1),RHOT CM EXPP=0.0 EXPP=1.0E-20 33 FORMAT(1X,'R,LAT,LNG,Z,T',1P5E12.3) 34 FORMAT(2X,'G1,RHOT',1P2E12.3) ENDIF !RT CM 2 DO 1 I=1,8 2 DO 1 I=1,IONS1 CC WRITE(6,1313) G(I),RHOT,EXPP,R(I) CC 1313 FORMAT(1X,'G,RHO,EX,R ',1P4E12.4) CM 1 R(I)=G(I)*RHOT*EXPP !ORIGINAL D&G VERSION RHOT=RHOT*CLTORFAC !MODIFIED JULY 26 1995 BY JDM 1 R(I)=G(I)*RHOT*EXPP C C FIX THE DENSITY DERIVATIVES C C DEFINE H AND Z WITH KM UNITS FOR DERIVATIVES C HM=H*PRAD ZM=Z*PRAD C CM DO JE=1,8 DO JE=1,IONS1 C DNDR=-2.0*RJ*((VLAT-.123*COSD(WLONG-21.0))/H)**2*R(JE) C DNDT=2.0*RJ*(RJ*VLAT-Z)*R(JE)/(H*H) C DNDPHI=2.0*RJ*RJ*(VLAT-(.123)*COSD(WLONG-21.0))*(.123)* DNDR=-2.0*Y1*((VLAT-.123*COSD(WLONG-21.0))/HM)**2*R(JE) DNDT=2.0*Y1*(Y1*VLAT-ZM)*R(JE)/(HM*HM) DNDPHI=2.0*Y1*Y1*(VLAT-(.123)*COSD(WLONG-21.0))*(.123)* #*SIND(WLONG-21.0)*R(JE)/(HM*HM) FP2(JE)=PLASCON(JE)*R(JE) DFP2DR(JE)=PLASCON(JE)*DNDR DFP2DT(JE)=PLASCON(JE)*DNDT DFP2DP(JE)=PLASCON(JE)*DNDPHI ENDDO !JE C C WRITE STATEMENTS C CX WRITE(6,*) ' DFPDR,DFPDT,DFPDP ',DFP2DR(1),DFP2DT(1),DFP2DP(1) C RETURN END SUBROUTINE INDSC(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON/DGPLAS/PLASCON,IONS1 DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8),PLASCON(8) DIMENSION R(8),G(8) DATA G/1.0,0.0,.07,.06,.06,.26,.06,.05/ RAD=.01745 PRAD=71300.0 Y1=RJ*PRAD H=1.82-.041*RJ Z=((7.*RJ-26.)/30.)*COS(RAD*(WLONG-21.)) RT=-((RJ*VLAT-Z)/H)**2 EXPP=0 IF(RT.LT.-80.)GO TO 2 EXPP=EXP(RT) 2 TEMP=100.-85.*EXPP V=(8.32*RJ+33.6) CM DO 1 I=1,8 DO 1 I=1,IONS1 1 R(I)=G(I)*RHOT*EXPP C C EVALUATE DENSITY DERIVATIVES C C C REDEFINE H AND Z WITH DIMENSIONS FOR DERIVATIVES C HM=H*PRAD ZM=Z*PRAD C C F1=RJ*VLAT-Z F1=Y1*VLAT-ZM CM DO JE=1,8 DO JE=1,IONS1 C DNDR=-2.0*(F1/H)*((VLAT-7.0*COSD(WLONG-21.0)/30.0)/ C #H+(F1*0.041)/(H*H))*R(JE) C DNDT=2.0*RJ*F1*R(JE)/(H*H) C DNDPHI=2.0*F1*(7.0*RJ-26.0)*SIND(WLONG-21.0)*R(JE)/ C #(30.0*H*H) DNDR=-2.0*(F1/HM)*((VLAT-7.0*COSD(WLONG-21.0)/30.0)/ #HM+(F1*0.041)/(HM*HM))*R(JE) DNDT=2.0*Y1*F1*R(JE)/(HM*HM) DNDPHI=2.0*F1*(7.0*Y1-26.0*PRAD)*SIND(WLONG-21.0)*R(JE)/ #(30.0*HM*HM) FP2(JE)=PLASCON(JE)*R(JE) DFP2DR(JE)=PLASCON(JE)*DNDR DFP2DT(JE)=PLASCON(JE)*DNDT DFP2DP(JE)=PLASCON(JE)*DNDPHI ENDDO !JE C C WRITE STATEMENTS C CX WRITE(6,*) ' DFPDR,DFPDT,DFPDP ',DFP2DR(1),DFP2DT(1),DFP2DP(1) C RETURN END SUBROUTINE WMTOR(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON/DGPLAS/PLASCON,IONS1 COMMON/BLKTOR/CLTORFAC,WMTORFAC,PLASFAC DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8),PLASCON(8) DIMENSION R(8),G(8) DATA G/1.0,0.0,.06,.08,.24,.25,.01,.01/ RAD=.01745 PRAD=71300.0 Y1=PRAD*RJ H=.2*SQRT(TEMP) V=12.6*RJ Z=RJ*.123*COS(RAD*(WLONG-21.)) EXPP=0 RT=-((RJ*VLAT-Z)/H)**2 IF(RT.LT.-80.) GO TO 2 EXPP=EXP(RT) CM 2 DO 1 I=1,8 2 DO 1 I=1,IONS1 CM 1 R(I)=G(I)*RHOT*EXPP !ORIGINAL D&G VERSION RHOT=RHOT*WMTORFAC !MODIFIED ON JULY 26, 1995 BY JDM 1 R(I)=G(I)*RHOT*EXPP C C EVALUATE THE DENSITY DERIVATIVES C C DEFINE H AND Z WITH KM UNITS FOR DERIVATIVES C ZM=Z*PRAD HM=H*PRAD C CM DO JE=1,8 DO JE=1,IONS1 C DNDR=-2.0*RJ*((VLAT-.123*COSD(WLONG-21.0))/H)**2*R(JE) C DNDT=2.0*RJ*(RJ*VLAT-Z)*R(JE)/(H*H) C DNDPHI=2.0*RJ*RJ*(VLAT-(.123)*COSD(WLONG-21.0))*(.123)* C #*SIND(WLONG-21.0)*R(JE)/(H*H) DNDR=-2.0*Y1*((VLAT-.123*COSD(WLONG-21.0))/HM)**2*R(JE) DNDT=2.0*Y1*(Y1*VLAT-ZM)*R(JE)/(HM*HM) DNDPHI=2.0*Y1*Y1*(VLAT-(.123)*COSD(WLONG-21.0))*(.123)* #*SIND(WLONG-21.0)*R(JE)/(HM*HM) FP2(JE)=PLASCON(JE)*R(JE) DFP2DR(JE)=PLASCON(JE)*DNDR DFP2DT(JE)=PLASCON(JE)*DNDT DFP2DP(JE)=PLASCON(JE)*DNDPHI ENDDO !JE C C WRITE STATEMENTS C CX WRITE(6,*) ' DFPDR,DFPDT,DFPDP ',DFP2DR(1),DFP2DT(1),DFP2DP(1) C RETURN END SUBROUTINE OTDSC(RJ,WLONG,VLAT,R,V,TEMP,RHOT, #DFP2DR,DFP2DT,DFP2DP,FP2) CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON/DGPLAS/PLASCON,IONS1 DIMENSION DFP2DR(8),DFP2DT(8),DFP2DP(8),FP2(8),PLASCON(8) DIMENSION R(8),G(8) DATA G/1.0,0.0,.07,.06,.06,.26,.06,.05/ RAD=.01745 PRAD=71300.0 Y1=PRAD*RJ H=1.0 Z=3.8*COS(RAD*(WLONG-21.-.9*(RJ-20.))) EXPP=0 RT=-((RJ*VLAT-Z)/H)**2 IF(RT.LT.-80.) GO TO 2 EXPP=EXP(RT) 2 TEMP=100.-85.*EXPP V=200. CM DO 1 I=1,8 DO 1 I=1,IONS1 1 R(I)=G(I)*RHOT*EXPP C C EVALUATE THE DENSITY DERIVATIVES C C C REDIFINE Z AND H WITH UNITS OF KM FOR DERIVATIVES C HM=H*PRAD ZM=Z*PRAD C F1=RJ*VLAT-Z F1=Y1*VLAT-ZM C F2=SIND((WLONG-21.0)-0.9*(RJ-20.0)) F2=SIND((WLONG-21.0)-0.9*(Y1-20.0*PRAD)/PRAD) CM DO JE=1,8 DO JE=1,IONS1 C DNDR=-2.0*(F1/(H*H))*(VLAT-20.0*0.19*0.9*F2)*R(JE) C DNDT=2.0*RJ*F1*R(JE)/(H*H) C DNDPHI=2.0*20.0*F1*0.19*F2*R(JE)/(H*H) DNDR=-2.0*(F1/(HM*HM))*(VLAT-20.0*0.19*0.9*F2)*R(JE) DNDT=2.0*Y1*F1*R(JE)/(HM*HM) DNDPHI=2.0*20.0*PRAD*F1*0.19*F2*R(JE)/(HM*HM) C WRITE(6,*) ' JE,PLASCON,R',JE,PLASCON(JE),R(JE) FP2(JE)=PLASCON(JE)*R(JE) DFP2DR(JE)=PLASCON(JE)*DNDR DFP2DT(JE)=PLASCON(JE)*DNDT DFP2DP(JE)=PLASCON(JE)*DNDPHI ENDDO !JE C C WRITE STATEMENTS C CX WRITE(6,*) ' DFPDR,DFPDT,DFPDP ',DFP2DR(1),DFP2DT(1),DFP2DP(1) C C C FOLLOWING LINES ARE FOR H+ C C EL=RJ/(COS(VLAT)**2) C Z0=20.*.123*COS(RAD*(WLONG-21.)-.9*(RJ-20.)) C ENP=(30.78-RJ)/16.9*(10.**-3) C R(2)=ENP*EXP(-ABS((RJ*VLAT-Z0)/2.)) C TYPE*,'L ',EL,'NP ',ENP RETURN END SUBROUTINE IONOS(HITE,THETA,ENION,FP2E,WIONOS,WMAGOS, #DWIDR,DWMDR) CM IMPLICIT REAL*16 (A-H,O-Z) C IMPLICIT REAL*8 (A-H,O-Z) C C C HASHIMOTO AND GOLDSTEIN IONOSPHERE DENSITY C C C THIS PROGRAM CALCULATES THE HASHIMOTO JOVIAN IONOSPHERE AND C DUMPS THE DATA TO SCREEN OR TO A FILE CALLED 'IONOS.DAT' C PI=3.1415927 DEGRAD=180.0/3.14159 THETA=THETA*DEGRAD COS2=COSD(THETA)*COSD(THETA) C DO J=1,NTIMES C HITE=2500.0-FLOAT(J-1)*50.0 C RO=71372.0/(SQRT(1.0-0.14371*COS2) RO=71300.0*(1.0-COS2/15.4) IF(HITE.GT.276.0) GO TO 2399 ENION=10.0 FP2E=8.0503D+07*ENION GO TO 7001 2399 VAL1=(991.0/(HITE-276.0))**4.31+1.0 PARM=1417.0-HITE APARM=ABS(PARM) CHECK=APARM/386.3 IF(CHECK.LT.293.0) VAL2=1.15E+05*EXP(PARM/386.3) IF(CHECK.GE.293.0) VAL2=0.0 ENION=VAL2/VAL1 !HASHIMOTO EXPRESSION FP2E=8.0503D+07*ENION 7001 CONTINUE C C DERIVATIVES OF THE DENSITY WRT R, THETA, AND PHI C C DEFINE CONSTANTS C AA=1.15D+05 A=1417.0 B=386.3 C=991.0 D=276.0 E=4.31 C FACT=1.0/(1.0+C/(HITE-D)**E) FACT2=E*C/((HITE-D)*(HITE-D)) FACT3=(C/(HITE-D))**E-1.0 FACT4=FACT*FACT DNIONDR=-AA*EXP((A-HITE)/B)*(FACT/B-FACT2*FACT3*FACT4) C C WE ASSUME DNIONDT=DNIONDP=0 C C DEFINE WEIGHTING FUNCTIONS C DEL=1.0 !CHANGES WEIGHTING FUNCTION WIDTH HC=2490.0-3.0*(90.0-THETA) !LINEAR FIT TO LATITUDE DEPENDENCE DH=HITE-HC HFUN=DH/DEL HDENOM=DEL*DEL+HFUN*HFUN WIONOS=0.5-ATAN(HFUN)/PI !IONOSPHERE WEIGHTING FUNCTION WMAGOS=0.5-ATAN(-HFUN)/PI !MAGNETOSPHERE WEIGHTING FUNC. DWIDR=-DEL/(PI*HDENOM) !DERIVATIVE OF IONOSPHERE FUNC. DWMDR=-DWIDR !DERIVATIVE OF MAGNETOSPHERE FUNC. CX WRITE(6,*) ' DH,DWIDR,WIONOS,WMAGOS ',DH,DWIDR,WIONOS,WMAGOS CX WRITE(6,*) ' H= ',HITE,' NIONOS= ',ENION RETURN END