C IRIFO7, DECEMBER 1979 C INTERNATIONAL REFERENCE IONOSPHERE (IRI). C THIS PROGRAM PRODUCES HEIGHT PROFILES OF TEMPERATURE C AND DENSITY OF THE ELECTRONS AND IONS IN THE IONOSPHERE FOR C SPECIFIED LOCATION,TIME AND SOLAR ACTIVITY. C ADDRESS: PROF.K.RAWER, HERRENSTR.43 C D 7801 MARCH , F.R.G. C INPUT: C 1.) LATITUDE(F6.1),LONGITUDE(F6.1),SOLAR-ACTIVITY(R)(F6.1), C MONTH(F4.1),HOUR(F4.1,1X),XHI(F6.1),BEGIN,END,STEPWIDTH(OF THE C REQUIRED HEIGHT RANGE)(3(F6.1,1X)),KOBE(I1,1X),JNE,JNEMAX, C JTN,JTE,JTI,JTETI,JO,JHHE,JO2,JNO,JMAG(11I1) C 2.) HMF2/KM AND NMF2/M-3 OR FOF2/MHZ (F5.1,E9.3) C IF HMF2 LESS 10.0 THEN THE CCIR VALUES FOR FOF2 AND HMF2 C ARE TAKEN BY USE OF THE CCIR-TAPE FOR FOF2 AND M3000. C THE USER SHOULD ADAPT PROCEDURE CCIRCA TO HIS OWN C CCIR-TAPE. C IF XHI LESS -10.0 THEN XHI WILL BE CALCULATED. C KOBE=0 SHOULD BE USED FOR PRINTER OUTPUT AND KOBE=1 FOR C PUNCHER OR TAPE OUTPUT.THE TEN INPUT VARIABLES C (JNE...JNO) ARE SWITCHED (1=YES,0=NO) TO CHOOSE YOUR C PARAMETERS.JMAG=0 MEANS GEOGRAFIC JMAG=1 GEOMAGNETIC C LATITUDE AND LONGITUDE. INTEGER EGNR,AGNR,SMONTH,DAYNR,DDO,IIF,DO2 REAL LATI,LONGI,MO2,MO,LOGE,MONTH,MODIP,NMF2,M3000,MAGBR, &NDEL0,NMF1,MUE,NME,NMD,K,NHABR,NDEL,NDX,NE0,NEI,MM,MLAT, &MLONG,NOBO2 DIMENSION F(3),B0F(2,2,8),OUTF(50,12),UF(38),RIF(4),XSM(2), &FF(38),FA(3,11),FZ(3,11),JF(10),DDO(4),IIF(4),DO2(2),MM(3), &Q(3),QT(3),PF1O(12),PF3O(12),HO(4),MO2(3),MO(5),G(144),E(4), &PF2O(4),HO2(2),CTNN(3),PG1O(80),PG2O(32),PG3O(80),CTN(3) LOGICAL WINTER,SUMMER,SCHALT,EXT,NIGHT,F1REG,VALLEY, &IRDLOW,IRDUPP,CCIREI,ANF COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR &/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2 COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR &/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2 &/BLOCK4/G/BLOCK5/ZX,TNX,ATN,CTN/BLOCK6/NIGHT,E &/BLOCK8/HS,TNHS,XSM,MM/BLO10/BETA,ETA,DELTA,ZETA EXTERNAL XE1,XE2,XE3,XE4,XE5,XE6,TEDER DATA FF,UF/4H(1X,,4HF6.1,1H,,35*1H ,4H(1H0,4H,3X,,4H1HH,, &3H2X,,34*1H /,FA/4H3X,E,4H10.4,1H,,4H3X,F,4H6.4,,1H , &4H3X,F,4H6.1,,1H ,4H3X,F,4H6.1,,1H ,4H3X,F,4H6.1,,1H , &4H3X,F,4H7.4,,1H ,4H3X,F,4H5.1,,1H ,4H3X,F,4H5.1,,1H , &4H3X,F,4H5.1,,1H ,4H3X,F,4H5.1,,1H ,4H3X,F,4H5.1,,1H /,FZ/ &4H7X,2,4HHNE,,3H3X,,4H4X,6,4HHN/N,4HMAX,,4H5X,2,4HHTN,,3H2X,, &4H5X,2,4HHTE,,3H2X,,4H5X,2,4HHTI,,3H2X,,4H4X,5,4HHTE/,3HTI,, &4H5X,4,4HHRDO,2H+,,4H4X,4,4HHRDH,2H+,,4H3X,5,4HHRDH,3HE+,, &4H3X,5,4HHRDO,3H2+,,4H3X,5,4HHRDN,3H0+,/,ENFF,ENUF/2*1H)/ LOGE=1/ALOG(10.0) UMR=0.01745329 PHI=3.1415927 CCIREI=.FALSE. C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS, C EGNR=INPUT, AGNR=OUTPUT, KONSOL=MONITOR OUTPUT EGNR=5 AGNR=6 KONSOL=6 READ(EGNR,2) LATI,LONGI,R,MONTH,HOUR,XHI,AH,EH,SH,KOBE,(JF(I), &I=1,10),JMAG,HMF2,FOF2 WRITE(KONSOL,2122) LATI,LONGI,R,MONTH,HOUR,XHI,AH,EH,SH, &KOBE,(JF(I),I=1,10),JMAG,HMF2,FOF2 2 FORMAT(3F6.1,2F4.1,F6.1,1X,3(F6.1,1X),I1,1X,11I1/F5.1,E9.3) 2122 FORMAT(1H*,3F6.1,2F4.1,F6.1,1X,3(F6.1,1X),I1,1X,11I1/F5.1,E9.3) IOND=JF(8)+JF(9)+JF(10)+JF(7) C C INPUT OF TRANSFORMATION COEFFICIENTS (PROCEDURE FIELDG) ...... C CALL KOEFFI(B0F) CALL KOEFFB(G) C C CALCULATION OF XHI,SUNSET,SUNRISE,COV,XHIM ................... C IF(JMAG.LE.0) GOTO 888 MLAT=LATI MLONG=LONGI 888 CALL GGM(JMAG,LONGI,LATI,MLONG,MLAT) DAYNR=IFIX(MONTH)*30-15 SUNDEC=-0.40915*COS(2.0*PHI/365.25*FLOAT(DAYNR+8)) Z1=SIN(SUNDEC)*SIN(LATI*UMR) Z2=COS(SUNDEC)*COS(LATI*UMR) IF(ABS(Z2).GT.0.0) GOTO 120 SAX=24.0 IF(Z1.GE.0.0) SAX=0.0 GOTO 140 120 IF(ABS(Z1/Z2).LE.1.0) GOTO 510 SAX=24.0 IF((Z1+Z2).GE.0.0) SAX=0.0 GOTO 140 510 SAX=12.0-ARCOS(-Z1/Z2)/(UMR*15.0) 140 SUX=24.0-SAX XLSTA=15.0*(HOUR-12.0) COV=63.75+R*(0.728+R*0.00089) IF(XHI.GT.(-10.0)) GOTO 520 COSXHI=Z1+Z2*COS(XLSTA*UMR) XHI=ARCOS(COSXHI)/UMR GOTO 503 520 COSXHI=COS(XHI*UMR) 503 XHI0=49.84733+0.349504*ABS(LATI) XHI100=38.96113+0.509932*ABS(LATI) XHIM=(XHI0+(XHI100-XHI0)*R/100.0)/UMR CALL FIELDG(LATI,LONGI,300.0,XMA,YMA,ZMA,BET,DIP,MODIP) MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR SUNDEC=SUNDEC/UMR C C CLASSIFICATION OF DIFFERENT TIMES AND REGIONS ................ C F1REG=.TRUE. VALLEY=.FALSE. SUMMER=.FALSE. WINTER=.FALSE. NIGHT=.FALSE. IF((MONTH.GT.10.0).OR.(MONTH.LT.3.0)) WINTER=.TRUE. IF((MONTH.GT.4.0).AND.(MONTH.LT.9.0)) SUMMER=.TRUE. EXT=SUMMER IF(LATI.GT.0.0) GOTO 1111 EXT=WINTER WINTER=SUMMER 1111 SUMMER=EXT IF((HOUR.GT.SUX).OR.(HOUR.LT.SAX)) NIGHT=.TRUE. IF((XHI.GT.XHIM).OR.NIGHT.OR.WINTER) F1REG=.FALSE. IF(NIGHT) F1REG=.FALSE. IF(WINTER) F1REG=.FALSE. C C INPUT OF THE ION DENSITY PARAMETER ARRAYS PF1O,PF2O AND PF3O...... C IF(IOND.LT.1) GOTO 141 IIF(1)=2 IF(ABS(LATI).LT.30.0) IIF(1)=1 IIF(2)=2 IF(COV.LT.100.0) IIF(2)=1 IIF(3)=3 IF(WINTER) IIF(3)=4 IF(SUMMER) IIF(3)=2 IIF(4)=1 IF(NIGHT) IIF(4)=2 CALL KOEFP1(PG1O) CALL KOEFP2 (PG2O) CALL KOEFP3 (PG3O) DO 6000 I=1, 4 RIF(I)=FLOAT(IIF(I)) 6000 CONTINUE CALL SUFE (PG1O,RIF,12,PF1O) CALL SUFE(PG2O,RIF,4,PF2O) CALL SUFE(PG3O,RIF,12,PF3O) C C CALCULATION OF ELECTRON DENSITY PARAMETERS................ C 141 DELL=4.32 IF(ABS(MODIP).GE.18.0) DELL=1.0+EXP(-(ABS(MODIP)-30.0)/10.0) IF(.NOT.F1REG) GOTO 150 C1=0.1244-(4.44E-4)*R+0.09/DELL FOF1=FOF1ED(MAGBR,R,COSXHI) NMF1=1.24E10*FOF1*FOF1 150 XHINON=XHIS(MONTH,12.0,LATI) FOE=FOEEDI(R,XHI,XHINON,LATI,HOUR,SUX,SAX) IF(HMF2.GE.10.0) GOTO 501 CALL F2OUT(MODIP,LATI,LONGI,MAGBR,R,MONTH,HOUR,FOE,HMF2,FOF2) GOTO 502 501 IF(FOF2.GT.100.0) FOF2=SQRT(FOF2/1.24E+10) 502 NMF2=1.24E10*FOF2*FOF2 NME=1.24E10*FOE*FOE NMD=XMDED(XHI,R,4.0E8) 7000 HME=HPOL(HOUR,110.0,105.0,SAX,SUX) HMD=HPOL(HOUR,81.0,88.0,SAX,SUX) SMONTH=IFIX(MONTH/3.0) IF(SMONTH.LT.1) SMONTH=4 B0=XPOL(B0F,HOUR,SAX,SUX,DELL,SMONTH,R) COS2=COS(MLAT*UMR) COS2=COS2*COS2 FLU=(COV-40.0)/30.0 ETA=0.058798-0.0070305*COS2+FLU*(-0.014065+0.0069724*COS2)+ &(0.0024287+0.0042810*COS2-0.00015280*FOF2)*FOF2 ZETA=0.078922-0.0046702*COS2+FLU*(-0.019132+0.0076545*COS2)+ &(0.0032513+0.0060290*COS2-0.00020872*FOF2)*FOF2 BETA=-128.03+20.253*COS2+FLU*(-8.0755-0.65896*COS2)+(0.44041 &+0.71458*COS2-0.042966*FOF2)*FOF2 XXX=EP1ST(-94.45,BETA) XXXX=DEP1ST(-94.45,BETA) DELTA=(ETA*XXX-ZETA/2.0)/(ETA*XXXX+ZETA/400.0) B1=3.0 F(1)=HPOL(HOUR,0.02+0.03/DELL,0.05,SAX,SUX) F(2)=HPOL(HOUR,4.6,4.5,SAX,SUX) F(3)=HPOL(HOUR,-11.5,-4.0,SAX,SUX) NDEL0=5.0 IF(WINTER) NDEL0=10.0 DNDH0=0.016 IF(SUMMER) DNDH0=0.01 NHABR=HPOL(HOUR,10.5/DELL,28.0,SAX,SUX) XXX=EPSTEP(45.0,67.0,-10.0,20.0,ABS(LATI)) HBR=HPOL(HOUR,17.8/DELL,XXX,SAX,SUX) NDEL=HPOL(HOUR,NDEL0/DELL,81.0,SAX,SUX) DNDHBR=HPOL(HOUR,DNDH0/DELL,0.06,SAX,SUX) IF(NDEL.GT.1.0) VALLEY=.TRUE. IF(.NOT.VALLEY) GOTO 600 CALL TAL(HME,NME,NHABR,NDEL,HBR,1.0,DNDHBR,EXT,E) IF(.NOT.EXT) GOTO 600 WRITE(KONSOL,650) 650 FORMAT(1X,"NO ELECTRON DENSITY E-VALLEY,AS THE & MODEL FUNCTION HAS A SECOND EXTREMUM IN THE VALLEY-REGION") VALLEY=.FALSE. 600 FP1=F(1) FP2=-FP1*FP1/2.0 FP30=(-F(2)*FP2-FP1+1.0/F(2))/(F(2)*F(2)) FP3U=(-F(3)*FP2-FP1-1.0/F(3))/(F(3)*F(3)) HDX=HMD+F(2) X=HDX-HMD NDX=NMD*EXP(X*(FP1+X*(FP2+X*FP30))) DNDX=NDX*(FP1+X*(2.0*FP2+X*3.0*FP30)) X=HME-HDX K=-DNDX*X/(NDX*ALOG(NDX/NME)) D1=DNDX/(NDX*K*X**(K-1.0)) IF(.NOT.VALLEY) HBR=0.0 IF(F1REG) GOTO 700 HMF1=0.0 NMF1=0.0 C1=0.0 700 HEF=HME+HBR C C SEARCH FOR HMF1,HST,HEF,AND HA.......................... C 7107 IF(.NOT.F1REG) NMF1=(NME+NMF2)/2.0 924 H=0.0 133 H=H+10.0 IF(H.GT.(HMF2-HME)) GOTO 135 CALL REGFA1(HMF2-10.0-H,HMF2,1.0,NMF1,XE2,SCHALT,HMF1) IF(SCHALT) GOTO 133 GOTO 137 135 WRITE(KONSOL,11) 11 FORMAT(1H ,1X," HMF1 IS & NOT EVALUATED BY THE FUNCTION XE2") HMF1=HMF2 C1=0.0 137 IF(HMF1.GT.(HEF+10.0)) GOTO 380 C1=0.0 HMF1=HMF2 WRITE(KONSOL,9) 9 FORMAT(1H ,1X," HEF GREATER AS HMF1") 380 H=0.0 125 H=H+3.0 IF(H.GT.(HMF1-HME)) GOTO 900 CALL REGFA1(HMF1,HME+H,1.0,NME,XE3,SCHALT,HST) STR=HST IF(SCHALT) GOTO 125 GOTO 360 900 IF (HMF2-HEF.GE.B0) GOTO 922 B0=(HMF2-HEF)/1.1 WRITE(KONSOL,923) B0 923 FORMAT(1X,"CORRECTION BY DECREASING OF B0 TO",F7.1) GOTO 924 922 B1=B1+0.5 IF(B1.GE.10.0) GOTO 901 H=0.0 GOTO 133 901 WRITE(KONSOL,100) 100 FORMAT(1H ,1X," HST IS NOT &EVALUATED BY THE FUNCTION XE3") HST=(HMF1+HME-HBR)/2.0 360 HZ=(HST+HMF1)/2.0 WRITE (KONSOL,902) B1 902 FORMAT (1X,"WE PUT B1= ", F3.1," TO GET HST") IF(HZ.GT.(HEF+10.0)) GOTO 950 HST=(HMF1+HEF)/2.0 HZ=(HST+HMF1)/2.0 950 D=HZ-HST HEF=HME+HBR T=D*D/(HZ-HEF-D) C C CALCULATION OF NEUTRAL TEMPERATURE PARAMETER........... C HTA=90.0 TNA=183.0 ZX=125.0 ANF=.TRUE. Z1=-180.0 882 TUN=TUNCAL(COV,LATI,SUNDEC,Z1) TNX=371.6678+0.0518806*TUN-294.3505*EXP(-0.00216222*TUN) ATN=0.63662*(TUN-TNX) HDEL=ZX-HTA TDEL=TNX-TNA HD2=HDEL*HDEL CTN(1)=1.9*TDEL/HDEL CTN(3)=3.0*TDEL/(HD2*HD2) CTN(2)=CTN(3)*1.333333*HDEL-CTN(1)/HD2 IF(.NOT.ANF) GOTO 881 ANF=.FALSE. TX=TN(130.0,TNX,ATN,CTN) TNXN=TNX ATNN=ATN CTNN(1)=CTN(1) CTNN(2)=CTN(2) CTNN(3)=CTN(3) Z1=XLSTA GOTO 882 C C CALCULATION OF ELECTRON TEMPERATURE PARAMETER......... C 881 HTA=120.0 H0T=200.0 H0N=400.0 IF(ABS(MLAT).GT.40.0) H0T=350.0 H0=H0T IF(NIGHT) H0=H0N HMAX=70.0*EXP((-1.4E-3)*MLAT*MLAT)+200.0 CALL TELAT(MLAT,1600.0,700.0,0.5,0.47,0.024,0.0,F1N,F2N) TX=TN(H0N,TNXN,ATNN,CTNN) TEN=TE(H0N,F1N,F2N,HMAX,0.0,H0N,0.0,0.0,0.0) DTEN=(TEN-TX)/(H0N-HTA) CALL TELAT(MLAT,2325.0,725.0,1.0,3.4,-0.014,2.56E4,F1T,F2T) TX=TN(H0,TNXN,ATNN,CTNN) XXX=TE(H0,F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0) XXXX=TE(H0,F1N,F2N,0.0,0.0,H0N,TX,HTA,DTEN) TE0=HPOL(HOUR,XXX,XXXX,SAX,SUX) TX=TN(H0,TNX,ATN,CTN) QUO=(TE0-TX)/(H0-HTA) C C CALCULATION OF ION TEMPERATURE PARAMETERS C 887 XSM(1)=430.0 MM(2)=HPOL(HOUR,3.0,0.0,SAX,SUX) Z1=EXP(-0.09*MLAT) YSM2=1240.0-1400.0*Z1/((1.0+Z1)*(1.0+Z1)) X2=ABS(MLAT) X1=X2*(0.47+X2*0.024)*UMR X3=COS(X1) X4=YSM2 TIN=1200.0-300.0*SIGN(1.0,X3)*SQRT(ABS(X3)) YSM2=TIN IF(X4.GE.TIN) YSM2=HPOL(HOUR,X4,TIN,SAX,SUX) TX=TN(XSM(1),TNXN,ATNN,CTNN) XXX=TE(XSM(1),F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0) XXXX=TE(XSM(1),F1N,F2N,0.0,0.0,H0N,TX,HTA,DTEN) Z1=HPOL(HOUR,XXX,XXXX,SAX,SUX) Z2=TN(XSM(1),TNX,ATN,CTN) IF(YSM2.LE.Z2) YSM2=(Z1+Z2)/2.0 CALL REGFA1(130.0,500.0,0.1,YSM2,TEDER,SCHALT,HS) IF(.NOT.SCHALT) GOTO 250 HS=200.0 250 TNHS=TN(HS,TNX,ATN,CTN) MM(1)=DTNDH(HS,ATN,CTN) MM(3)=MM(2) XSM(2)=XSM(1) C C SMOOTHING OF THE ION TEMPERATURE TO KEEP IT LESS THAN C THE ELECTRON TEMPERATURE AT HSM C HSM=1000.0 XXX=TE(HSM,F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0) XXXX=TE(HSM,F1N,F2N,HMAX,0.0,H0N,0.0,0.0,0.0) TESM=HPOL(HOUR,XXX,XXXX,SAX,SUX) IF(TESM.GE.TI(HSM)) GOTO 240 XXX=DTEDH(HSM,F2T,HMAX,2.0,H0T,0.0,0.0) XXXX=DTEDH(HSM,F2N,HMAX,0.0,H0N,0.0,0.0) MM(3)=HPOL(HOUR,XXX,XXXX,SAX,SUX) XSM(2)=(TESM-10.0-YSM2+MM(2)*XSM(1)-MM(3)*HSM)/(MM(2)-MM(3)) C C CALCULATION OF ION DENSITY PARAMETER.................. C 240 IF(IOND.LT.1) GOTO 189 Z1=0.0 IF(XHI.LE.90.0) Z1=COSXHI HFIXO=300.0 IF((IIF(2).EQ.2).AND.(IIF(3).EQ.2)) HFIXO=249.0 MO(1)=EPSTEP(PF1O(1),PF1O(2),PF1O(3),PF1O(4),Z1) MO(2)=EPSTEP(PF1O(5),PF1O(6),PF1O(7),PF1O(8),Z1) MO(3)=0.0 HO(1)=EPSTEP(PF1O(9),PF1O(10),PF1O(11),PF1O(12),Z1) HO(4)=PF2O(1) MO(4)=PF2O(2) MO(5)=PF2O(3) DDO(1)=9 DDO(2)=5 DDO(3)=5 DDO(4)=50 7100 HO(2)=290.0 IF((IIF(2).EQ.2).AND.(IIF(3).EQ.2)) HO(2)=237.0 HO(3)=(4.60517-MO(5)*(HO(4)-PF2O(4)))/MO(4)+HO(4) IF(HO(2).LT.HO(3)) GOTO 7101 MO(4)=MO(4)-0.001 GOTO 7100 7101 Z2=5.0 X=HO(2) Z1=0.0 7102 X=X+Z2 Y=RPID(X,HFIXO,98.0,4,MO,DDO,HO) IF(Y.LE.Z1) GOTO 7103 Z1=Y GOTO 7102 7103 IF(Z2.LE.1.0) GOTO 7104 X=X-Z2 Z2=1.0 GOTO 7102 7104 H0O=X-0.5 IF(XHI.GE.90.0) COSXHI=0.0 DO 7105 I=2,4,2 L=I/2 HO2(L)=PF3O(1+I)+PF3O(2+I)*COSXHI 7105 MO2(L+1)=PF3O(7+I)+PF3O(8+I)*COSXHI DO2(1)=5 DO2(2)=5 MO2(1)=PF3O(7)+PF3O(8)*COSXHI 7106 Y=RPID(H0O,PF3O(1),PF3O(2),2,MO2,DO2,HO2) IF(Y.LE.0.1) GOTO 189 MO2(3)=MO2(3)-0.02 GOTO 7106 C C CALCULATION FOR THE REQUIRED HEIGHT RANGE....................... C 189 HA=65.0 IF(NIGHT) HA=80.0 IF(AH.LT.HA) AH=HA IF(EH.GT.1000.0) EH=1000.0 KOMB=JF(4)+JF(5)+JF(6) IF(((EH-AH)/SH+1.0).LT.50.0) GOTO 230 EH=AH+49.0*SH WRITE(AGNR,190) EH 190 FORMAT(1H ,1X," TOO MANY HEIGHT STEPS,EH IS REDUCED & TO ",F4.0) 230 I=0 X=EH 300 I=I+1 IN=1 OUTF(I,IN)=X IF((JF(1)+JF(2)).LE.0) GOTO 330 NEI=XE(X) IF(JF(1).LE.0) GOTO 340 IN=IN+1 OUTF(I,IN)=NEI 340 IF(JF(2).LE.0) GOTO 330 IN=IN+1 OUTF(I,IN)=NEI/NMF2 330 IF(KOMB.LE.0) GOTO 7108 X1=TN(X,TNX,ATN,CTN) IF(X.GE.HTA) GOTO 7109 Z2=-1.0 TEH=-1.0 GOTO 7110 7109 TX=TN(X,TNXN,ATNN,CTNN) Z1=TE(X,F1N,F2N,HMAX,0.0,H0N,TX,HTA,DTEN) Z2=X1 IF(X.GE.HS) Z2=TI(X) IF(X.LT.H0) GOTO 7119 XXX=TE(X,F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0) TEH=HPOL(HOUR,XXX,Z1,SAX,SUX) GOTO 7110 7119 TEH=TE(X,0.0,0.0,0.0,0.0,H0,X1,HTA,QUO) 7110 IF(JF(3).LE.0) GOTO 7112 IN=IN+1 OUTF(I,IN)=X1 7112 IF(JF(4).LE.0) GOTO 7113 IN=IN+1 OUTF(I,IN)=TEH 7113 IF(JF(5).LE.0) GOTO 7120 IN=IN+1 OUTF(I,IN)=Z2 7120 IF(JF(6).LE.0) GOTO 7108 IN=IN+1 OUTF(I,IN)=SIGN(1.0,TEH)*TEH/Z2 7108 IF(IOND.LE.0) GOTO 7118 Z1=RPID(X,HFIXO,98.0,4,MO,DDO,HO) Z2=RPID(X,PF3O(1),PF3O(2),2,MO2,DO2,HO2) C FIRST COMPUTE NOBO2,THE RATIO OF NO+ TO O2+ AT HEIGHT OF REL. C O+ MAX.(H0O).SAME RATIO IS CONSERVED AT ALL HEIGHTS ABOVE IF(JF(7).LE.0) GOTO 7114 Z1B= RPID (H0O,HFIXO,98.0,4,MO,DDO,HO) Z2B= RPID (H0O,PF3O(1),PF3O(2),2,MO2,DO2,HO2) NOBO2= (100.0-Z1B-Z2B)/Z2B IF(JF(7).LE.0) GOTO 7114 IN=IN+1 OUTF(I,IN)=Z1 7114 IF(JF(8).LE.0) GOTO 7116 IN=IN+1 CALL RDHHE(X,H0O,Z1,Z2,NOBO2,10.0,OUTF(I,IN),OUTF(I,IN+1)) 7116 IF(JF(9).LE.0) GOTO 7117 IN=IN+2 OUTF(I,IN)=Z2 7117 IF(JF(10).LE.0) GOTO 7118 IN=IN+1 OUTF(I,IN)=RDNO(X,H0O,Z2,Z1,NOBO2) 7118 X=X-SH IF(X.GE.AH) GOTO 300 IEI=I C C OUTPUT ON THE SPECIFIED DEVICE................................. C IF(KOBE.NE.1) GOTO 7020 WRITE(AGNR,7030) (JF(I),I=1,10),IEI,LATI,LONGI,R,MONTH,HOUR 7030 FORMAT(1X,10I1,2X,I3,2(2X,F6.1),2X,F5.1,2X,F4.1,2X,F5.2) 7020 IF(KOBE.NE.0) GOTO 7041 WRITE (AGNR,7051) 7051 FORMAT (1X,"INPUT@") IF (JMAG.LT.1) GOTO 7053 WRITE (AGNR,7052) MLAT,MLONG 7052 FORMAT (1X,"MLAT=",F6.1,2X,"MLONG=",F6.1) GOTO 7055 7053 WRITE (AGNR,7054) LATI,LONGI 7054 FORMAT (1X,"LATI=",F6.1,2X,"LONGI=",F6.1) 7055 WRITE(AGNR,7050) R,MONTH,HOUR 7050 FORMAT(1X,"R=",F5.0," MONTH=",F4.1," HOUR=",F5.2) WRITE (AGNR,7061) 7061 FORMAT (1X,"CALCULATED VALUES@") IF (JMAG.LT.1) GOTO 7063 WRITE (AGNR,7062) LATI,LONGI 7062 FORMAT (1X,"LATI=",F6.1,2X,"LONGI=",F6.1) GOTO 7065 7063 WRITE (AGNR,7064)MLAT,MLONG 7064 FORMAT (1X,"MLAT=",F6.1,2X,"MLONG=",F6.1) 7065 WRITE(AGNR,7060) DIP,MODIP,MAGBR,XHI 7060 FORMAT(1X,"DIP=",F6.1," MODIP=",F6.1," MAGLA=",F6.1, &" XHI=",F6.1) WRITE(AGNR,7066) SAX,SUX,SUNDEC 7066 FORMAT(1X,"SUNRISE",F4.1,"L.T. SUNSET@",F4.1,"L.T.", &5X,"SUNDEC.=",F6.1) Z1=0.0 IF(F1REG) Z1=HMF1 Z2=0.0 IF(F1REG) Z2=NMF1 WRITE(AGNR,7070) NMF2,Z2,NME,NMD 7070 FORMAT(1X,"NMF2=",E7.2," NMF1= ",E7.2," NME=",E7.2, &" NMD=",E7.2) WRITE(AGNR,7080) HMF2,Z1,HME,HMD 7080 FORMAT(1X,"HMF2=",F5.1," HMF1=",F5.1," HME=",F5.1," HMD=",F5.1) KARL=4 DO 7091 J=1,7 IF(JF(J).LT.1) GOTO 7091 DO 7092 I=1,3 FF(KARL+I)=FA(I,J) 7092 UF(KARL+I)=FZ(I,J) KARL=KARL+3 7091 CONTINUE J=8 IF(JF(J).LT.1) GOTO 7093 DO 7094 I=1,3 FF(KARL+I)=FA(I,J) UF(KARL+I)=FZ(I,J) II=I+3 JJ=J+1 FF(KARL+II)=FA(I,JJ) 7094 UF(KARL+II)=FZ(I,JJ) KARL=KARL+6 7093 CONTINUE DO 7095 J=9, 10 IF(JF(J).LT.1) GOTO 7095 JJ=J+1 DO 7096 I=1,3 FF(KARL+I)=FA(I,JJ) 7096 UF(KARL+I)=FZ(I,JJ) KARL=KARL+3 7095 CONTINUE FF(KARL+1)=ENFF UF(KARL+1)=ENUF WRITE(AGNR,UF) 7041 DO 7040 I=1,IEI 7040 WRITE(AGNR,FF) (OUTF(I,L),L=1,IN) IF(KOBE.EQ.1) WRITE(AGNR,35) 35 FORMAT(1X,"111111") STOP END C D.BILITZA,H.THIEMANN,K.RAWER,IPW FREIBURG, C AUG.79 ********************** C IRI-PROCEDURES *************************************** C \REF.@ K. RAWER, S. RAMAKRISHNAN, D. BILITZA, C INTERNATIONAL REFERENCE IONOSPHERE 1978, U.R.S.I. BRUSSELS' C FOR CALCULATION OF ELECTRON DENSITY, TEMPERATURE, ION TEMPERATURE C AND ION RELATIVE DENSITY IN THE HEIGHT RANGE 70 - 1000 KM.? C I. ELECTRON DENSITY ---------------------------------- FUNCTION EP1ST(X,BET) C EPSTEIN-STEP-FUNCTION WITH TRANSITION FROM Y=0 TO Y=1 AT X=0 C AND TRANSITION THICKNESS BET EP1ST=1.0/(1.0+EXP(-X/BET)) RETURN END FUNCTION DEP1ST(X,BET) U=EXP(-X/BET) DEP1ST=U/(BET*(1.0+U)*(1.0+U)) RETURN END FUNCTION XE1(H) C REPRESENTING ELECTRON DENSITY C FOR HEIGHTS NOT GREATER 1000 KM AND NOT LESS HMF2 C BY HARMONIZED BENT-MODEL ADMITTING VARIABILITY OF C GLOBAL PARAMETERS@ ETA, ZETA, BETA, DELTA WITH GEOM. LATITUDE C (MLAT), SOLAR FLUX(FLU) AND CRITICAL FREQUENCY(FOF2). C ALO GLOBAL ARE PEAK-PARAMETERS NMF2,HMF2? C \REF.@K.RAWER,S.RAMAKRISHNAN,1978' REAL NMF2 COMMON/BLOCK1/HMF2,NMF2,HMF1/BLO10/BETA,ETA,DELTA,ZETA X=(H-HMF2)/(1000.0-HMF2)*700.0+300.0-DELTA Y=(1000.0-HMF2)/700.0*(BETA*ETA*ALOG((1.0+EXP((X-394.50)/ &BETA))/(1.0+EXP((-94.5-DELTA)/BETA)))+ZETA*(100.0*ALOG((1.0+ &EXP((X-300.0)/100.0))/(1.0+EXP(-DELTA/100.0)))-X+300.0-DELTA)) XE1=NMF2*EXP(-Y) RETURN END REAL FUNCTION XE2(H) C ELECTRON DENSITY FOR HEIGHTS LESS HMF2 AND NOT LESS HMF1 REAL NMF2 COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR X=(HMF2-H)/B0 IF(ABS(X).LT.1.0E-10) GOTO 100 XE2=NMF2*EXP(-X**B1)/COSH(X) 100 XE2=NMF2 200 RETURN END REAL FUNCTION XE3(H) C ELECTRON DENSITY FOR HEIGHTS LESS HMF1 AND NOT LESS HZ REAL NMF2 COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR XE3=XE2(H)+NMF2*C1*SQRT(ABS(HMF1-H)/B0) RETURN END REAL FUNCTION XE4(H) C ELECTRON DENSITY FOR HEIGHTS LESS HZ AND NOT LESS HEF COMMON/BLOCK2/B0,B1,C1,HZ,T,HST,STR/BLOCK3/ &HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2 A=((STR/HST-1)*(H-HZ))/(HEF-HZ)+1 XE4=XE3 (A*(HZ+T/2.0-SIGN(1.0,T)*SQRT(T*(HZ-H+T/4.0)))) RETURN END REAL FUNCTION XE5(H) C ELECTRON DENSITY FOR HEIGHTS LESS HEF AND NOT LESS HME(VALLEY) REAL NME,NMD,K LOGICAL NIGHT COMMON/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1, &FP2/BLOCK6/NIGHT,E(4) T3=H-HME T1=T3*T3*(E(1)+T3*(E(2)+T3*(E(3)+T3*E(4)))) IF(NIGHT) GOTO 100 XE5=NME*(T1+1.0) GOTO 200 100 XE5=NME*EXP(T1) 200 RETURN END REAL FUNCTION XE6(H) C ELECTRON DENSITY FOR HEIGHTS LESS HME AND NOT LESS HA REAL NME,NMD,K COMMON/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2 IF(H.GT.HDX) GOTO 100 Z=H-HMD FP3=FP3U IF(Z.GT.0.0) FP3=FP30 XE6=NMD*EXP(Z*(FP1+Z*(FP2+Z*FP3))) GOTO 200 100 Z=HME-H XE6=NME*EXP(-D1*Z**K) 200 RETURN END REAL FUNCTION XE(H) C ELECTRON DENSITY BEETWEEN HA(KM) AND 1000 KM C SUMMARIZING PROCEDURES NE1....6? REAL NMF2,NME,NMD,K COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR &/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2 IF(H.LT.HMF2) GOTO 100 XE=XE1(H) GOTO 200 100 IF(H.LT.HMF1) GOTO 300 XE=XE2(H) GOTO 200 300 IF(H.LT.HZ) GOTO 400 XE=XE3(H) GOTO 200 400 IF(H.LT.HEF) GOTO 500 XE=XE4(H) GOTO 200 500 IF(H.LT.HME) GOTO 600 XE=XE5(H) GOTO 200 600 XE=XE6(H) 200 RETURN END C II. ELECTRON TEMPERATURE ---------------------------------- FUNCTION DTEDH(H,F2,HMAX,D,H0,DTNDH,A) C D.BILITZA 1978, CALCULATES THE HEIGHT DERIVATIVE C ELECTRON TEMPERATURE PROFILE AFTER PROCEDURE TE IF(H.LT.H0) GOTO 100 Y=EXP(-0.03*(H-HMAX)) F=1.0+Y DTEDH=D-F2*0.03*Y*(1.0-Y)/(F*F*F) RETURN 100 DTEDH=DTNDH+A RETURN END SUBROUTINE TELAT(PHI,A,B,XN,P1,P2,C,F1,F2) C D. BILITZA 1978, CALCULATES THE PARAMETERS F1 AND F2 C NEEDED IN THE ELECTRON TEMPERATURE MODEL PROCEDURE TE, C BOTH DEPENDING ON GEOMAGNETIC LATITUDE (PHI) XPHI=ABS(PHI) F=(P1+P2*XPHI)*XPHI F=F*0.01745329 Y=COS(F) F1=A-B*SIGN(1.0,Y)*ABS(Y)**XN Y=EXP(-0.1*XPHI) F=1.0+Y F2=C*Y/(F*F) RETURN END FUNCTION TE(H,F1,F2,HMAX,D,HH0,TNH,HTA,A) C D. BILITZA 1978. ELECTRON TEMPERATURE PROFILE C IN THE HEIGHT RANGE 120-1000 KM. F1,F2 ARE LATITUDE DEPENDENT AND C ARE GIVEN BY TELAT. HMAX IS THE HEIGHT OF A POSSIBLE MAXIMUM C D THE TEMPERATURE HEIGHT GRADIENT, H0 THE INTERSECTION HEIGHT C BELOW WHICH THE PROFILE IS BUILT TO APPROACH THE CIRA 72 TEMP. C TNA AT HTA . A TNH ARE COEFFICIENTS OF THE FITTING FUNCTION IF(H.LT.HH0) GOTO 10 Y=EXP(-0.03*(H-HMAX)) F=1.0+Y TE=F1+F2*Y/(F*F)+D*(H-700.0) RETURN 10 TE=TNH+A*(H-HTA) RETURN END C III. ION TEMPERATURE -------------------------------------- FUNCTION TUNCAL(COV,XLATI,SD,SLSTA) C CALCULATES THE EXOSPHERIC TEMPERATURE FOR COVINGTON-INDEX C COV,LATITUDE XLATI,SOLARDEKLINATION SD AND LOCAL SOLAR TIME C ANGLE SLSTA USING CIRA72-MODEL UMR=0.01745329 TC=379.0+3.24*COV ETA=ABS(XLATI-SD)/2.0 THETA=ABS(XLATI+SD)/2.0 H1=COS(ETA*UMR)**2.2 H2=SIN(THETA*UMR)**2.2 TD=1.0+0.3*H1 TN=1.0+0.3*H2 A=(TD-TN)/TN TAU=SLSTA-37.0+6.0*SIN((SLSTA+43.0)*UMR) X=COS(TAU/2.0*UMR) TUNCAL=TC*TN*(1.0+A*SIGN(1.0,X)*ABS(X)**3.0) RETURN END REAL FUNCTION TN(X,TTNX,ATN,CCTN) C NEUTRAL TEMPERATURE C D. BILITZA, 1978, NEUTRAL TEMPERATURE PROFILE C APPROXIMATING ROUGHLY CIRA72-MODEL? DIMENSION CCTN(3) Z=X-125.0 IF(Z.LE.0.0) GOTO 100 Y=Z**2.5 Y=(1.0+(4.5E-6)*Y)*Z YY=CCTN(1)/ATN Y=YY*Y Y=ATAN(Y)*ATN TN=TTNX+Y GOTO 200 100 TN=TTNX+Z*(CCTN(1)+Z*Z*(CCTN(2)+Z*CCTN(3))) 200 RETURN END FUNCTION DTNDH(H,ATN,CCTN) C DERIVATIVE OF NEUTRAL TEMPERATURE C D.BILITZA, 1978. APPROXIMATE HEIGHT DERIVATIVE OF NEUTRAL C CIRA-72 MODEL TEMPERATURE PROFILE DIMENSION CCTN(3) Z=H-125.0 IF(Z.GT.0.0) GOTO 100 DTNDH=CCTN(1)+Z*Z*(3.0*CCTN(2)+Z*4.0*CCTN(3)) RETURN 100 H1=CCTN(1)/ATN H2=Z**2.5 H3=H1*Z*(1.0+(4.5E-6)*H2) DTNDH=ATN/(1.0+H3*H3)*H1*(1.0+(15.75E-6)*H2) RETURN END REAL FUNCTION TI(H) C ION TEMPERATURE FOR HEIGHTS NOT GREATER 1000 KM AND NOT LESS HS REAL MM,G(2) COMMON/BLOCK8/HS,TNHS,XSM(2),MM(3) G(1)=20.0 G(2)=50.0 SUM=MM(1)*(H-HS)+TNHS DO 100 I=1,2 A=H-XSM(I) IF((H-XSM(I)).LT.(100.0*G(I))) A=G(I)*ALOG(1.0+ &EXP((H-XSM(I))/G(I))) B=HS-XSM(I) IF((HS-XSM(I)).LT.(100.0*G(I))) B=G(I)*ALOG(1.0+ &EXP((HS-XSM(I))/G(I))) 100 SUM=SUM+(MM(I+1)-MM(I))*(A-B) TI=SUM RETURN END REAL FUNCTION TEDER(H) C THIS FUNCTION ALONG WITH PROCEDURE REGFA1 ALLOWS TO FIND C THE HEIGHT ABOVE WHICH TN BEGINS TO BE DIFFERENT FROM TI REAL MM COMMON/BLOCK5/ZX,TNX,ATN,CTN(3)/BLOCK8/HS,TNHS,XSM(2),MM(3) TNH=TN(H,TNX,ATN,CTN) DTDX=DTNDH(H,ATN,CTN) TEDER=DTDX*(XSM(1)-H)+TNH RETURN END C IV. ION RELATIVE PRECENTAGE DENSITY REAL FUNCTION RPID (H, H0, N0, M, ST, ID, XS) C D.BILITZA,1977,THIS ANALYTIC FUNCTION IS USED TO REPRESENT THE C RELATIVE PRECENTAGE DENSITY OF ATOMAR AND MOLECULAR OXYGEN IONS. C THE M+1 HEIGHT GRADIENTS ST(M+1) ARE CONNECTED WITH EPSTEIN- C STEP-FUNCTIONS AT THE STEP HEIGHTS XS(M) WITH TRANSITION C THICKNESSES ID(M). RPID(H0,H0,N0,....)=N0. C INSTEAD OF 88.0 YOU CAN USE THE HIGHEST ALLOWED ARGUMENT C FOR EXP AT YOUR COMPUTER. REAL N0 DIMENSION ID(4), ST(5), XS(4) SUM=(H-H0)*ST(1) DO 100 I=1,M A=H-XS(I) XI=FLOAT(ID(I)) IF (A.LT.88.0*XI) A=XI*ALOG(1.0+EXP((H-XS(I))/XI)) B=H0-XS(I) IF (B.LT.88.0*XI) B=XI*ALOG(1.0+EXP((H0-XS(I))/XI)) 100 SUM=SUM+(ST(I+1)-ST(I))*(A-B) SUM=N0*EXP(SUM) IF (SUM.LT.0.00005) SUM=0.0 RPID=SUM RETURN END SUBROUTINE RDHHE (H,HB,RDOH,RDO2H,RNO,PEHE,RDH,RDHE) C RAWER, OCT.79,H+ AND HE+ RELATIVE PERECENTAGEDENSITY FOR HEIGHTS C NOT GREATER THAN 1000 KM.THE O+ AND O2+ REL. PER. DENSITIES SHOULD C BE GIVEN (RDOH,RDO2H). HB IS THE HEIGHT WERE N(O+)=100_. C RNO IS THE RATIO OF NO+ TO O2+DENSITY AT H=HB. IF (H-HB)100,100,200 100 RDH=0.0 GOTO 300 200 RDH=(100.0-RDOH-(1.0+RNO)*RDO2H)*(1.0-PEHE/100.0) 300 RDHE=RDH*PEHE/(100.0-PEHE) RETURN END REAL FUNCTION RDNO(H,HB,RDO2H,RDOH,RNO) C D.BILITZA, 1978. NO+ RELATIVE PERCENTAGE DENSITY FOR HEIGHTS C NOT LESS THAN 100 KM.FOR MORE INFORMATION SEE SUBROUTINE RDHHE. IF (H-HB) 100,100,200 100 Y=100.0-RDO2H-RDOH GOTO 300 200 Y=RNO*RDO2H 300 IF(Y.LT.0.0005) Y=0.0 RDNO=Y RETURN END C END IRI-FUNCTIONS****************************************** REAL FUNCTION XMDED (XHI,R,XW) C D. BILITZA, 1978, CALCULATES ELECTRON DENSITY OF D MAXIMUM. C XHI IS SOLAR ZENITH ANGLE, R ZURICH SUNSPOT NUMBER AND XW THE C REQUIRED NIGHT VALUE \SOURCE@MECHTLY-BILITZA,1974'. Y=6.05+0.088*R YW=XW/1.0E8 Z=(-0.1/(ALOG(YW/Y)))**0.3704 SUXHI=ARCOS(Z) IF (SUXHI.LT.1.0472) SUXHI=1.0472 XXHI=XHI/57.2957795 IF (XXHI.GT.SUXHI) GOTO 100 X=COS(XXHI) XMDED=Y*1.0E8*EXP(-0.1/X**2.7) RETURN 100 XMDED=YW*1.0E8 RETURN END REAL FUNCTION FOEEDI(R,XHI,XHIM,SLATI,ST,SU,SA) C D.BILITZA, 1977, CALCULATES FOE BY THE EDINBURGH-METHOD. C INPUT@ZUERICH SUNSPOT NUMBER(R),SOLAR ZENITH ANGLE AT LOCAL C HOUR ST (XHI) AND AT NOON (XHIM), LATITUDE (SLATI),SUNSET C AND SUNRISE (SU, SA/HOUR). C \REF.@ KOURIS-MUGGELETON, CCIR DOC. 6/3/07 SEPT.73' XLATI=ABS(SLATI) COV=63.75+R*(0.728+0.00089*R) A=1.0+0.0094*(COV-66.0) SL=COS(XLATI*0.0174533) IF(XLATI.LT.32.0) GOTO 100 SM=0.11-0.49*SL C=92.0+35.0*SL GOTO 400 100 SM=-1.93+1.92*SL C=23.0+116.0*SL 400 IF(XHIM.GT.89.0) XHIM=89.0 B=COS(XHIM*0.0174533)**SM SP=1.31 IF(XLATI.GT.12.0) SP=1.2 DX=0.0 IF(XHI.GT.73.0.AND.XHI.LE.90.0) DX=(6.27E-13)*(XHI-50.0)**8.0 IF(ST.GT.SU) GOTO 300 IF(ST.GT.SA) GOTO 600 D=0.077**SP*EXP(-1.68*(SA-ST)) GOTO 200 600 D=COS((XHI-DX)*0.0174533)**SP GOTO 200 300 D=0.077**SP*EXP(-1.01*(ST-SU)) 200 FOE=A*B*C*D SP=1.0+0.0098*R SMIN=0.017*SP*SP IF(FOE.LT.SMIN) FOE=SMIN FOEEDI=(FOE)**0.25 RETURN END REAL FUNCTION FOF1ED(YLATI,R,COSXHI) C D. BILITZA,1977, CALCULATES FOF1 FOR DIP- LATITUDE YLATI , C ZURICH SUNSPOT NUMBER R AND COSINUS OF SOLAR ZENITH C ANGLE COSXHI. \REF@E.D.DUCHARME ET.AL.,RADIO SCI.,8,10,837-839, C HOWEVER WITH MAG. INSTEAD OF GEOM. LATITUDE,R.EYFRIG,1979' XLATI=ABS(YLATI) F0=4.35+XLATI*(0.0058-(1.2E-4)*XLATI) F100=5.348+XLATI*(0.011-(2.3E-4)*XLATI) FS=F0+(F100-F0)*R/100.0 XMUE=0.093+XLATI*(0.0046-(5.4E-5)*XLATI)+(3.0E-4)*R FOF1ED=FS*COSXHI**XMUE RETURN END REAL FUNCTION HMF2ED(XMAGBR,R,X,XM3) C D. BILITZA,CALCULATES THE PEAK HEIGHT HMF2 FOR THE MAGNETIC C LATITUDE XMAGBR AND THE ZUERICH SUNSPOT NUMBER R C USING M3000(XM3) AND THE RATIO FOF2/FOE. C \REF. D. BILITZA ET. AL. TELECOMM.J.,46,549-553,1979' F1=(2.32E-3)*R+0.222 F2=1.2-(1.16E-2)*EXP((2.39E-2)*R) F3=0.096*(R-25.0)/150.0 DELM=F1*(1.0-R/150.0*EXP(-XMAGBR*XMAGBR/1600.0))/(X-F2)+F3 HMF2ED=1490.0/(XM3+DELM)-176.0 RETURN END SUBROUTINE REGFA1(X11,X22,EPS,FW,F,SCHALT,X) C REGULA-FALSI-PROCEDURE TO FIND X WITH F(X)-FW=0. X1,X2 ARE THE C STARTING VALUES. THE COMUTATION ENDS WHEN THE X-INTERVAL C HAS BECOME LESS THAN EPS . IF SIGN(F(X1)-FW)= SIGN(F(X2)-FW) C THEN SCHALT=.TRUE. LOGICAL L1,LINKS,K,SCHALT SCHALT=.FALSE. X1=X11 X2=X22 F1=F(X1)-FW F2=F(X2)-FW IF(ABS(SIGN(1.0,F1)-SIGN(1.0,F2))-10E-10) 100,110,110 110 K=.FALSE. NG=2 200 X=(X1*F2-X2*F1)/(F2-F1) GOTO 400 300 L1=LINKS IF(LINKS) GOTO 500 X=X1+FLOAT(NG-1)*(X2-X1)/FLOAT(NG) GOTO 400 500 X=X1+(X2-X1)/FLOAT(NG) 400 FX=F(X)-FW LINKS=SIGN(1.0,F1)*SIGN(1.0,FX).GT.0 K=.NOT.K IF(LINKS) GOTO 600 X2=X F2=FX GOTO 700 600 X1=X F1=FX 700 IF(ABS(X2-X1).LE.EPS) GOTO 800 IF(K) GOTO 300 IF((LINKS.AND.(.NOT.L1)).OR.(.NOT.LINKS.AND.L1)) NG=2*NG GOTO 200 100 X=0.0 SCHALT=.TRUE. 800 RETURN END SUBROUTINE TAL(SHMAX,STMAX,SHABR,SDELTA,SHBR,STE0,SDTDH0,AUS6,SPT) C D. BILITZA,15.1.77,CAL. OF THE COEF. SPT(4) FOR A POLYNOM OF C 5TH DEGREE (SPT(0)=1) TO FIT THE VALLEY REPR. BY THE PARAMETER@ C SHMAX IS THE HEIGHT OF THE LOWER MAXIMUM WHERE THE VALLEY STARTS, C STMAX IS THE DENSITY AT SHMAX, SHABR IS THE HEIGHT DIFFERENCE C BETWEEN THE HEIGHT OF THE DEEPEST VALLEY DENSITY AND SHMAX, SDELTA C IS THE PERCENTAGE DEPTH,SHBR THE WIDTH,STE0 THE UPPER FIXPOINT C NORMALISED TO STMAX AND SDTDH0 IS THE LOGARITHMIC DERIV. AT C THIS FIXPOINT.IF THERE IS A MINIMUM OR A SECOND MAXIMUM IN THE C VALLEY REGION AUS6=.TRUE. DIMENSION SPT(4) LOGICAL AUS6 X=SHABR Y=SHBR IF(SDELTA.LE.60.0) GOTO 100 Z1=ALOG(1.0-SDELTA/100.0)/(X*X) Z2=ALOG(STE0)/(Y*Y) Z3=SDTDH0/(2.0*Y) GOTO 200 100 Z1=-SDELTA/(X*X*100.0) Z2=(STE0-1.0)/(Y*Y) Z3=SDTDH0*STE0/(2.0*Y) 200 Z4=X-Y SPT(4)=2.0*(Z1*(Y-2.0*X)*Y-Z2*(X-2.0*Y)*X+Z3*Z4*X)/(X*Y*Z4*Z4*Z4) SPT(3)=(Z1*(2.0*Y-3.0*X)+Z2*X)/(X*Z4*Z4)-(2.0*X+Y)*SPT(4) SPT(2)=-2.0*Z1/X-2.0*X*SPT(3)-3.0*X*X*SPT(4) SPT(1)=Z1-X*(SPT(2)+X*(SPT(3)+X*SPT(4))) AUS6=.FALSE. Z1=(4.0*SPT(3)+5.0*SPT(4)*X)/(10.0*SPT(4)) Z2=Z1*Z1+2.0*SPT(1)/(5.0*SPT(4)*X) IF(Z2.LT.0.0) GOTO 300 Z3=SQRT(Z2) Z2=-1.0*Z1+Z3 IF(Z2.GT.0.0.AND.Z2.LT.Y) AUS6=.TRUE. Z2=-1.0*Z1-Z3 IF(Z2.GT.0.0.AND.Z2.LT.Y) AUS6=.TRUE. 300 RETURN END REAL FUNCTION XHIS(XMONTH,XHOUR,XLATI) C CALCULATES THE SOLAR-ZENITH-ANGLE DAYNR=XMONTH*30.0-15.0 SUNDEC=-0.40915*COS(1.7202E-2*(DAYNR+8.0)) UMR=1.7453E-2 COSXHI=SIN(XLATI*UMR)*SIN(SUNDEC) +COS(XLATI*UMR)*COS(SUNDEC)* &COS(2.6180E-1*(XHOUR-12.0)) XHIS=ARCOS(COSXHI)/UMR RETURN END SUBROUTINE GGM(ART,XLG,BG,XLM,BM) C CALCULATES GEOM. COORDINATES (XLM,BM) FROM GEOG. COORDINATES C (XLG,BG) FOR ART=0 AND REVERSE FOR ART=1 INTEGER ART ZPI=6.2831853 FAKTOR=0.0174533 CBG=11.4*FAKTOR CI=COS(CBG) SI=SIN(CBG) IF(ART.EQ.0) GOTO 10 CBM=COS(BM*FAKTOR) SBM=SIN(BM*FAKTOR) CLM=COS(XLM*FAKTOR) SLM=SIN(XLM*FAKTOR) SBG=SBM*CI-CBM*CLM*SI BG=ARSIN(SBG) CBG=COS(BG) SLG=(CBM*SLM)/CBG CLG=(SBM*SI+CBM*CLM*CI)/CBG XLG=ARCOS(CLG) IF(SLG.LT.0.0) XLG=ZPI-ARCOS(CLG) BG=BG/FAKTOR XLG=XLG/FAKTOR XLG=XLG-69.8 IF(XLG.LT.0.0) XLG=XLG+360.0 GOTO 20 10 YLG=XLG+69.8 CBG=COS(BG*FAKTOR) SBG=SIN(BG*FAKTOR) CLG=COS(YLG*FAKTOR) SLG=SIN(YLG*FAKTOR) SBM=SBG*CI+CBG*CLG*SI BM=ARSIN(SBM) CBM=COS(BM) SLM=(CBG*SLG)/CBM CLM=(-SBG*SI+CBG*CLG*CI)/CBM XLM=ARCOS(CLM) IF(SLM.LT.0.0) XLM=ZPI-ARCOS(CLM) BM=BM/FAKTOR XLM=XLM/FAKTOR 20 RETURN END SUBROUTINE KOEFFB(FIELD) C SHEIK,1977,TRANSFORMATION COEFF. G(144) VALID FOR 1973 FOR C CALCULATING THE MAGNETIC FIELD ACCORDING TO POGO 68/10 DIMENSION FIELD (144) REAL FEL1 (72),FEL2(72) DATA FEL1/0.0, 0.1506723,0.0101742, -0.0286519, 0.0092606, & -0.0130846, 0.0089594, -0.0136808,-0.0001508, -0.0093977, & 0.0130650, 0.0020520, -0.0121956, -0.0023451, -0.0208555, & 0.0068416,-0.0142659, -0.0093322, -0.0021364, -0.0078910, & 0.0045586, 0.0128904, -0.0002951, -0.0237245,0.0289493, & 0.0074605, -0.0105741, -0.0005116, -0.0105732, -0.0058542, &0.0033268, 0.0078164,0.0211234, 0.0099309, 0.0362792, &-0.0201070,-0.0046350,-0.0058722,0.0011147,-0.0013949, & -0.0108838, 0.0322263, -0.0147390, 0.0031247, 0.0111986, & -0.0109394,0.0058112, 0.2739046, -0.0155682, -0.0253272, & 0.0163782, 0.0205730, 0.0022081, 0.0112749,-0.0098427, & 0.0072705, 0.0195189, -0.0081132, -0.0071889, -0.0579970, & -0.0856642, 0.1884260,-0.7391512, 0.1210288, -0.0241888, & -0.0052464, -0.0096312, -0.0044834, 0.0201764, 0.0258343, &0.0083033, 0.0077187/ DATA FEL2/0.0586055,0.0102236,-0.0396107, & -0.0167860, -0.2019911, -0.5810815,0.0379916, 3.7508268, & 1.8133030, -0.0564250, -0.0557352, 0.1335347, -0.0142641, & -0.1024618,0.0970994, -0.0751830,-0.1274948, 0.0402073, & 0.0386290, 0.1883088, 0.1838960, -0.7848989,0.7591817, & -0.9302389,-0.8560960, 0.6633250, -4.6363869, -13.2599277, & 0.1002136, 0.0855714,-0.0991981, -0.0765378,-0.0455264, & 0.1169326, -0.2604067, 0.1800076, -0.2223685, -0.6347679, &0.5334222, -0.3459502,-0.1573697, 0.8589464, 1.7815990, &-6.3347645, -3.1513653, -9.9927750,13.3327637, -35.4897308, &37.3466339, -0.5257398, 0.0571474, -0.5421217, 0.2404770, & -0.1747774,-0.3433644, 0.4829708,0.3935944, 0.4885033, & 0.8488121, -0.7640999, -1.8884945, 3.2930784,-7.3497229, & 0.1672821,-0.2306652, 10.5782146, 12.6031065, 8.6579742, & 215.5209961, -27.1419220,22.3405762,1108.6394043/ K=0 DO 10 I=1,72 K=K+1 FIELD(K)=FEL1(I) 10 FIELD(72+K)=FEL2(I) RETURN END SUBROUTINE FIELDG(DLAT,DLONG,ALT,X,Y,Z,F,DIP,SMODIP) C THIS IS A SPECIAL VERSION OF THE POGO68/107 C MAGNETIC FIELD LEGENDRE MODEL (IN UNITS OF GAUSS). C F IS TOTAL FIELD, Z THE (DOWNWARD) VERTICAL COMPONENT WHILE C X AND Y ARE COMPONENTS IN THE EQUAT. PLANE X TO ZERO LONGITUDE). C DIP=INCLINATION ANGLE/DEGREE. MODIP=RAWER"S MODFIED DIP. C REFERENCE@ METEOROLOGICAL AND ASTRONOMICAL INFLUENCES ON RADIO C WAVE PROPAGATION. PERGAMON PRESS,1963,ALSO RADIO SCI. 10,669,75 C INPUT@ DLAT, DLONG=GEOGRAPHIC COORDINATES/DEG., ALT=ALTITUDE/KM. C THE SET OF COEFFICIENTS IS GIVEN BY SUBROUTINE KOEFFB. DIMENSION H(144),XI(3) COMMON/BLOCK4/G(144) RLAT=DLAT*0.0174533 CT=SIN(RLAT) ST=COS(RLAT) NMAX=11 D=SQRT(40680925.0-272336.0*CT*CT) RLONG=DLONG*0.0174533 CP=COS(RLONG) SP=SIN(RLONG) ZZZ=(ALT+40408589.0/D)*CT/6371.2 RHO=(ALT+40680925.0/D)*ST/6371.2 XXX=RHO*CP YYY=RHO*SP RQ=1.0/(XXX*XXX+YYY*YYY+ZZZ*ZZZ) XI(1)=XXX*RQ XI(2)=YYY*RQ XI(3)=ZZZ*RQ IHMAX=NMAX*NMAX+1 LAST=IHMAX+NMAX+NMAX IMAX=NMAX+NMAX-1 DO 100 I=IHMAX,LAST 100 H(I)=G(I) DO 200 K=1,3,2 I=IMAX IH=IHMAX 300 IL=IH-I F1=2.0/FLOAT(I-K+2) X1=XI(1)*F1 Y1=XI(2)*F1 Z1=XI(3)*(F1+F1) I=I-2 IF((I-1).LT.0) GOTO 400 IF((I-1).EQ.0) GOTO 500 DO 600 M=3,I,2 H(IL+M+1)=G(IL+M+1)+Z1*H(IH+M+1)+X1*(H(IH+M+3)-H(IH+M-1))- &Y1*(H(IH+M+2)+H(IH+M-2)) H(IL+M)=G(IL+M)+Z1*H(IH+M)+X1*(H(IH+M+2)-H(IH+M-2))+ &Y1*(H(IH+M+3)+H(IH+M-1)) 600 CONTINUE 500 H(IL+2)=G(IL+2)+Z1*H(IH+2)+X1*H(IH+4)-Y1*(H(IH+3)+H(IH)) H(IL+1)=G(IL+1)+Z1*H(IH+1)+Y1*H(IH+4)+X1*(H(IH+3)-H(IH)) 400 H(IL)=G(IL)+Z1*H(IH)+2.0*(X1*H(IH+1)+Y1*H(IH+2)) 700 IH=IL IF(I.GE.K) GOTO 300 200 CONTINUE S=0.5*H(1)+2.0*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2)) XT=(RQ+RQ)*SQRT(RQ) X=XT*(H(3)-S*XXX) Y=XT*(H(4)-S*YYY) Z=XT*(H(2)-S*ZZZ) F=SQRT(X*X+Y*Y+Z*Z) Z=-Z*CT-(Y*SP+X*CP)*ST DIP=ARSIN(Z/F) SMODIP=ARSIN(DIP/SQRT(DIP*DIP+ST)) DIP=DIP*57.2957795 9 SMODIP=SMODIP*57.2957795 RETURN END SUBROUTINE F2OUT(XMODIP,XLATI,XLONGI,XMAGBR,R,XMONTH,HOUR,FOE, &XHMF2,FOF2) C D.BILITZA,20.8.78,CALCULATES XHMF2 AND NMF2 USING THE C CCIR-MAPS FOR FOF2 AND M3000.YOU HAVE TO ADJUST THE TAPE C READING PROCEDURE CCIRCA TO YOUR COMPUTER SYSTEM AND CCIR-TAPE INTEGER QF(9),QM(7) DIMENSION FF0(988) REAL MM0(441) DATA QF/11,11,8,4,1,0,0,0,0/,QM/6,7,5,2,1,0,0/ REWIND 5 LMONTH=IFIX(XMONTH) CALL CCIRCA(R,LMONTH,FF0,MM0) FOF2=GAMMA1(XMODIP,XLATI,XLONGI,.FALSE.,HOUR,6,QF,9,76,13,988, &FF0) XM3000=GAMMA1(XMODIP,XLATI,XLONGI,.FALSE.,HOUR,4,QM,7,49,9,441, &MM0) XHMF2=HMF2ED(XMAGBR,R,FOF2/FOE,XM3000) RETURN END REAL FUNCTION GAMMA1(SMODIP,SLAT,SLONG,UT,XHOUR,IHARM,NQ,K1,M, &MM,M3,SFE) C SHEIKH,4.3.77,CALCULATES THE RATIO (FOF2/M3000) USING CCIR C NUMERICAL MAP COEFFICIENTS SFE(M,2*IHARM+1). C NQ(K1) IS AN INTEGER ARRAY GIVING NQ1,NQ2,NQ3,...AS HIGHEST C DEGREE OF LATITUDE FOR EACH LONGITUDE HARMONIC (IHARM). C M=1+NQ1+2(NQ2+1)+2(NQ3+1)+... .IF UT IS TRUE,THEN HOUR IS C TAKEN AS UNIVERSAL TIME,OTHERWISE IT IS ASSUMED AS LOCAL C TIME. DIMENSION NQ(K1),XSINX(13),COEF(100),C(12),S(12),SFE(M3) LOGICAL UT RD=57.2957795 IF(UT) GOTO 100 TIME=(15.0*XHOUR-180.0-SLONG)/RD GOTO 150 100 TIME=(15.0*XHOUR-180.0)/RD 150 S(1)=SIN(TIME) C(1)=COS(TIME) DO 250 I=2,IHARM C(I)=C(1)*C(I-1)-S(1)*S(I-1) S(I)=C(1)*S(I-1)+S(1)*C(I-1) 250 CONTINUE DO 300 I=1,M COEF(I)=SFE((I-1)*MM+1) DO 300 J=1,IHARM COEF(I)=COEF(I)+SFE((I-1)*MM+2*J)*S(J)+SFE((I-1)*MM+2*J+1)*C(J) 300 CONTINUE SUM=COEF(1) SS=SIN(SMODIP/RD) S3=SS XSINX(1)=1.0 INDEX=NQ(1) DO 350 J=1,INDEX SUM=SUM+COEF(1+J)*SS XSINX(J+1)=SS SS=SS*S3 350 CONTINUE XSINX(NQ(1)+2)=SS*S3 NP=NQ(1)+1 SS=COS(SLAT/RD) S3=SS DO 400 J=2,K1 S0=SLONG*FLOAT(J-1)/RD S1=COS(S0) S2=SIN(S0) INDEX=NQ(J)+1 DO 450 L=1,INDEX NP=NP+1 SUM=SUM+COEF(NP)*XSINX(L)*SS*S1 NP=NP+1 SUM=SUM+COEF(NP)*XSINX(L)*SS*S2 450 CONTINUE SS=SS*S3 400 CONTINUE GAMMA1=SUM RETURN END SUBROUTINE CCIRCA(R,IMONTH,FOF2,SM3000) C THIS SUBROUTINE CALCULATES THE COEFFICIENTS ARRAYS FOF2(76,13) AND C SM3000(49,9) FOR SOLARACTIVITY (R) AND MONTH USING OUR SPECIAL C CCIR-FORTRAN-TAPE. DIMENSION FOF2(988),SM3000(441),F2(13,76,2),FM3(9,49,2) DO 10 I=1,IMONTH READ(10,1) F2 10 READ(10,1) FM3 1 FORMAT(8E15.8) DO 20 I=1,76 DO 20 J=1,13 K=J+13*(I-1) 20 FOF2(K)=(F2(J,I,1)*(100.0-R)+F2(J,I,2)*R)/100.0 DO 30 I=1,49 DO 30 J=1,9 K=J+9*(I-1) 30 SM3000(K)=(FM3(J,I,1)*(100.0-R)+FM3(J,I,2)*R)/100.0 RETURN END SUBROUTINE KOEFFI(BOF) C DATA FOR THICKNESS OF BOTTOMSIDE F2 REGION ELECTRON DENSITY. DIMENSION BOF(2,2,8) REAL FELD(32) DATA FELD/114.0,64.0,134.0,77.0,128.0,66.0,75.0,73.0,113.0,115.0, &150.0,116.0,138.0,123.0,94.0,132.0,72.0,84.0,83.0,89.0,75.0,85.0, &57.0,76.0,102.0,100.0,120.0,110.0,107.0,103.0,76.0,86.0/ L=0 K=0 DO 10 I=1,2 DO 10 J=1,2 DO 10 K=1,8 L=L+1 10 BOF(I,J,K)=FELD(L) RETURN END SUBROUTINE KOEFP1(PG1O) C THIEMANN,1979,COEFFICIENTS PG1O FOR CALCULATING O+ PROFILES C BELOW THE F2-MAXIMUM. CHOSEN TO APPROACH DANILOV- C SEMENOV"S COMPILATION. DIMENSION PG1O(80) REAL FELD (80) DATA FELD/-11.0,-11.0,4.0,-11.0,0.08018, &0.13027,0.04216,0.25 ,-0.00686,0.00999, &5.113,0.1 ,170.0,180.0,0.1175,0.15,-11.0, &1.0 ,2.0,-11.0,0.069,0.161,0.254,0.18,0.0161, &0.0216,0.03014,0.1,152.0,167.0,0.04916, &0.17,-11.0,2.0,2.0,-11.0,0.072,0.092,0.014,0.21, &0.01389,0.03863,0.05762,0.12,165.0,168.0,0.008, &0.258,-11.0,1.0,3.0,-11.0,0.091,0.088, &0.008,0.34,0.0067,0.0195,0.04,0.1,158.0,172.0, &0.01,0.24,-11.0,2.0,3.0, -11.0,0.083,0.102, &0.045,0.03,0.00127,0.01,0.05,0.09,167.0,185.0, &0.015,0.18/ K=0 DO 10 I=1,80 K=K+1 10 PG1O(K)=FELD(I) RETURN END SUBROUTINE KOEFP2(PG2O) C THIEMANN,1979,COEFFICIENTS FOR CALCULATION OF O+ PROFILES C ABOVE THE F2-MAXIMUM (DUMBS,SPENNER@AEROS-COMPILATION) DIMENSION PG2O(32) REAL FELD(32) DATA FELD/1.0,-11.0,-11.0,1.0,695.0,-.000781, &-.00264,2177.0,1.0,-11.0,-11.0,2.0,570.0, &-.002,-.0052,1040.0,2.0,-11.0,-11.0,1.0,695.0, &-.000786,-.00165,3367.0,2.0,-11.0,-11.0,2.0, &575.0,-.00126,-.00524,1380.0/ K=0 DO 10 I=1,32 K=K+1 10 PG2O(K)=FELD(I) RETURN END SUBROUTINE KOEFP3(PG3O) C THIEMANN,1979,COEFFICIENTS FOR CALCULATING O2+ PROFILES. C CHOSEN AS TO APPROACH DANILOV-SEMENOV"S COMPILATION. DIMENSION PG3O(80) REAL FELD(80) DATA FELD/-11.0,1.0,2.0,-11.0,160.0,31.0,130.0, &-10.0,198.0,0.0,0.05922,-0.07983, &-0.00397,0.00085,-0.00313,0.0,-11.0,2.0,2.0,-11.0, &140.0,30.0,130.0,-10.0, &190.0,0.0,0.05107,-0.07964,0.00097,-0.01118,-0.02614, &-0.09537, &-11.0,1.0,3.0,-11.0,140.0,37.0,125.0,0.0,182.0, &0.0,0.0307,-0.04968,-0.00248, &-0.02451,-0.00313,0.0,-11.0,2.0,3.0,-11.0, &140.0,37.0,125.0,0.0,170.0,0.0, &0.02806,-0.04716,0.00066,-0.02763,-0.02247,-0.01919, &-11.0,-11.0,4.0,-11.0,140.0,45.0,136.0,-9.0, &181.0,-26.0,0.02994,-0.04879, &-0.01396,0.00089,-0.09929,0.05589/ K=0 DO 10 I=1,80 K=K+1 10 PG3O(K)=FELD(I) RETURN END SUBROUTINE SUFE (FIELD,RFE,M,FE) C SPECIAL SUBR. FOR CHOOSING THE REQUIRED ION DENSITY PARAMETER C ARRAY. THE INPUT FIELD INCLUDES DIFF. SETS OF DIMENSION M EACH C CARACTERISED BY 4 HEADER NUMBERS. RFE(4) SHOULD CONTAIN THE C CHOSEN HEADER NUMBERS.FE(M) IS THE CORRESPONDING SET. DIMENSION RFE(4),FE(12),FIELD(80),EFE(4) K=0 100 DO 101 I=1,4 K=K+1 101 EFE(I)=FIELD(K) DO 111 I=1,M K=K+1 111 FE(I)=FIELD(K) DO 120 I=1,4 IF((EFE(I).GT.-10.0).AND.(RFE(I).NE.EFE(I))) GOTO 100 120 CONTINUE RETURN END FUNCTION HPOL(XHOUR,TW,XNW,SA,SU) C D.BILITZA,1978 PROCEDURE FOR TIME-INTERPOLATION C USING EPSTEIN STEP FUNCTION OF ONE HOUR WIDTH (DAY/NIGHT C FUNCTION). TW,NW ARE THE DAY AND NIGHT VALUE OF THE PARAMETER. C SA, SU ARE TIME (HOURS) OF SUNRISE AND SUNSET/HOUR? IF(ABS(SU-SA).GE.24.0) GOTO 100 HPOL=XNW+(TW-XNW)/(1.0+EXP(-(XHOUR-SA)/1.0))+(XNW-TW) &/(1.0+EXP(-(XHOUR-SU)/1.0)) GOTO 200 100 HPOL=XNW IF(SA.LT.1.0) HPOL=TW 200 RETURN END REAL FUNCTION XPOL(AA,XHOUR,SA,SU,DELL,IMONTH,R) C SPECIAL FUNCTION TO INTERPOLATE ARRAY AA(LATI,R,MONTH,HOUR)= C A(2,2,4,2) IN LATITUDE (LATI), SOLAR ACTIVITY (R) AND TIME (HOUR) DIMENSION AA(2,2,8),SIPH(2),SIPL(2) JMONTH=IMONTH*2 DO 10 ISR=1,2 DO 15 ISLAT=1,2 15 SIPH(ISLAT)=HPOL(XHOUR,AA(ISLAT,ISR,JMONTH-1), &AA(ISLAT,ISR,JMONTH),SA,SU) SIPL(ISR)=SIPH(1)+(SIPH(2)-SIPH(1))/DELL 10 CONTINUE XPOL=SIPL(1)+(SIPL(2)-SIPL(1))/90.0*(R-10.0) RETURN END FUNCTION EPSTEP(BE,AB,TH,FIX,COSI) C EPSTEIN STEP FUNCTION = TRANSITION WITH THICKNESS TH BETWEEN C BE AND AB CENTERED WHERE VARIABLE COSI EQUALS FIX. EPSTEP=AB+(BE-AB)/(1.0+EXP(-(COSI-FIX)/TH)) RETURN END