C FILE USER1:[FACULTY.LEVIS]GROUNDE.FOR;1 5/14/90 c C FIELD INTENSITY CALCULATION AFTER BREMMER-VAN DER POL. C AFTER H.BREMMER,TERRESTRIAL RADIO WAVES,PP.104-110 AND 124 C AND K.A.NORTON,THE CALCULATION OF GROUND-WAVE FIELD INTENSITY OVER C A FINITELY CONDUCTING SPHERICAL EARTH,PROC.IRE,1941,P.623-639 C AND F.E.TERMAN,RADIO ENGINEERS HANDBOOK,PP.676,677 C PROGRAM BY DR MIROSLAV JOACHIM,C.C.I.R.,GENEVA REAL K,LAMBDA COMPLEX DELTA,J COMMON/SOL/AR,ARL,D,DA,DB,DELTA,EPS,F,J,K,LAMBDA,PI,PSI,SIG,WRK 1,ARHO AR=10.0 ARL=1.0 PI=3.14159265358979 J=(0.,1.) write(*,*)'Ground wave program for Pt=1 kW Vertical Polarization' write(*,*)'Enter real dielectric constant, conductivity (s/m),' write(*,*)'and frequency (MHz) :' read(*,*) EPS,SIEM,F write(*,*)'Output file ground.out contains R in km,' write(*,*)'spherical distance in km, field intensity in mV/m,' write(*,*)'and dB relative to one microvolt/m' open(unit=1,file='ground.out',status='unknown') SIG=1.E-11*SIEM LAMBDA=300./F WRK=6.E12*SIG*LAMBDA PSI=ATAN(EPS/WRK)-0.5*ATAN((EPS-1.)/WRK) K=0.002924*(LAMBDA**(1./3.))*SQRT(EPS*EPS+WRK*WRK)/ *SQRT(SQRT((EPS-1.)**2.+WRK*WRK)) DELTA=K*CEXP(J*(0.75*PI-PSI)) A=5.*LAMBDA**(1./3.) B=10.*LAMBDA**(1./3.) DO 19 ID=1,10 D=ID IF(D-B) 17,19,19 17 CONTINUE CALL FIELD1(E1) WRITE(1,18)DA,D,E1,DB 18 FORMAT(2F10.1,2F20.8) 19 CONTINUE DO 22 ID=2,10 DI=ID D=10.*DI IF(D-B) 21,22,22 21 CONTINUE CALL FIELD1(E1) WRITE(1,18)DA,D,E1,DB 22 CONTINUE DO 24 ID=2,19 DI=ID D=100.*DI IF(D-B) 23,24,24 23 CONTINUE CALL FIELD1(E1) WRITE(1,18)DA,D,E1,DB 24 CONTINUE DO 27 ID=1,10 DI=ID D=10.*DI IF(D-A) 27,27,26 26 CONTINUE CALL FIELD2(E2) WRITE(1,18)DA,D,E2,DB 27 CONTINUE DO 29 ID=2,10 DI=ID D=100.*DI IF(D-A) 29,29,28 28 CONTINUE CALL FIELD2(E2) IF(E2.LT.1.E-06) GO TO 32 WRITE(1,18)DA,D,E2,DB 29 CONTINUE DO 34 ID=12,18,2 DI=ID D=100.*DI IF(D-A) 34,34,35 35 CONTINUE CALL FIELD2(E2) IF(E2.LT.1.E-06) GO TO 32 WRITE(1,18)DA,D,E2,DB 34 CONTINUE DO 31 ID=2,10 DI=ID D=1000.*DI IF(D-A) 31,31,30 30 CONTINUE CALL FIELD2(E2) IF(E2.LT.1.E-06) GO TO 32 WRITE(1,18)DA,D,E2,DB 31 CONTINUE 32 STOP END C------------------------------------------------ SUBROUTINE FIELD1(E1) REAL K,LAMBDA COMPLEX Y1 COMPLEX DELTA,E,J,RHO,Y COMMON/SOL/AR,ARL,D,DA,DB,DELTA,EPS,F,J,K,LAMBDA,PI,PSI,SIG,WRK 1,ARHO DCIR=D DARC=12756.*SIN(DCIR/12756.) D=DARC DA=DCIR RHO=3140.*D*SQRT((EPS-1.)**2.+WRK*WRK)*CEXP(2.*J*PSI)/((EPS*EPS+WR 1K*WRK)*LAMBDA) ARHO=CABS(RHO) IF(ARHO-AR) 2,1,1 1 CONTINUE Y=-1./2./RHO Y1=Y DO 6 I=2,100 AI=2*(I-1)+1 Y1=Y1*AI/2./RHO AY1=CABS(Y1) IF(AY1.LT.1.E-04) GO TO 11 6 Y=Y+Y1 11 CONTINUE GO TO 3 2 CONTINUE IF(ARHO-ARL) 4,5,5 4 CONTINUE X=WRK B2=ATAN(EPS/X) B1=ATAN((EPS-1.)/X) B=2.*B2-B1 P=1000.*PI*D*COS(B2)*COS(B2)/X/LAMBDA/COS(B1) A=(2.+0.3*P)/(2.+P+0.6*P*P)-SQRT(P/2.)*EXP(-5.*P/8.)*SIN(B) E1=300.*A/D GO TO 7 5 CONTINUE Y1=2.*RHO Y=1.+J*CSQRT(PI*RHO)*CEXP(-RHO)-Y1 DO 10 I=2,200 AI=2*(I-1)+1 Y1=Y1*2.*RHO/AI AY1=CABS(Y1) IF(AY1.LT.1.E-04) GO TO 3 Y=Y+(-1)**I*Y1 10 CONTINUE 3 CONTINUE E=300.*(Y+CSQRT(RHO)*(-J*0.886*DELTA**3+CEXP(J*PI/6.)*1.62*DELTA** 14-CEXP(-J*PI/6.)*DELTA**5*(0.455+0.304*RHO)))/D E1=CABS(E) 7 CONTINUE DB=20.*ALOG10(1000.*E1) RETURN END c------------------------------------------------ SUBROUTINE FIELD2(E2) DIMENSION TAUI(15),TAUR(15) REAL K,LAMBDA COMPLEX DELTA,E(15),J,TAU(15) COMMON/SOL/AR,ARL,D,DA,DB,DELTA,EPS,F,J,K,LAMBDA,PI,PSI,SIG,WRK 1,ARHO DCIR=D DARC=12756.*SIN(DCIR/12756.) D=DARC DA=DCIR WRK2=0.25*PI+PSI WRK3=75.*PI/180.+3.*PSI WRK4=75.*PI/180.-5.*PSI WRK5=15.*PI/180.-PSI WRK6=60.*PI/180.-4.*PSI WRK7=4.*PSI WRK8=2.*PSI WRK9=15.*PI/180.+3.*PSI CHI=0.0537*D/(LAMBDA**(1./3.)) IF(K-0.6) 1,4,4 1 CONTINUE TAUR(1)=0.928+K*COS(WRK2)+1.237*K**3.*COS(WRK3)-0.5*K**4.*COS(WRK7 1)-2.755*K**5.*COS(WRK4) TAUR(2)=1.622+K*COS(WRK2)+2.163*K**3.*COS(WRK3)-0.5*K**4.*COS(WRK7 1)-8.422*K**5.*COS(WRK4) TAUR(3)=2.191+K*COS(WRK2)+2.921*K**3.*COS(WRK3)-0.5*K**4.*COS(WRK7 1)-15.36*K**5.*COS(WRK4) TAUI(1)=1.607-K*SIN(WRK2)-1.237*K**3.*SIN(WRK3)+0.5*K**4.*SIN(WRK7 1)-2.755*K**5.*SIN(WRK4) TAUI(2)=2.810-K*SIN(WRK2)-2.163*K**3.*SIN(WRK3)+0.5*K**4.*SIN(WRK7 1)-8.422*K**5.*SIN(WRK4) TAUI(3)=3.795-K*SIN(WRK2)-2.921*K**3.*SIN(WRK3)+0.5*K**4.*SIN(WRK7 1)-15.36*K**5.*SIN(WRK4) DO 2 I=4,15 AI=I-1 TAUR(I)=1.116*(AI+0.75)**(2./3.)+K*COS(WRK2) TAUI(I)=1.932*(AI+0.75)**(2./3.)-K*SIN(WRK2) 2 CONTINUE GO TO 3 4 CONTINUE TAUR(1)=0.4043+0.618*COS(WRK5)/K-0.236*SIN(WRK8)/K**2.-0.0533*COS( 1WRK9)/K**3.+0.00226*COS(WRK6)/K**4. TAUR(2)=1.288+0.1940*COS(WRK5)/K-0.0073*SIN(WRK8)/K**2.+0.0120*COS 1(WRK9)/K**3.-0.00160*COS(WRK6)/K**4. TAUI(1)=0.7003-0.6183*SIN(WRK5)/K+0.2364*COS(WRK8)/K**2.-0.0533*SI 1N(WRK9)/K**3.-0.00226*SIN(WRK6)/K**4. TAUI(2)=2.232-0.1940*SIN(WRK5)/K+0.0073*COS(WRK8)/K**2.+0.0120*SIN 1(WRK9)/K**3.+0.00160*SIN(WRK6)/K**4. DO 5 I=3,15 AI=I-1 TAUR(I)=1.116*(AI+0.25)**(2./3.)+0.2241*COS(WRK5)/(K*(AI+0.25)**(2 1./3.)) TAUI(I)=1.932*(AI+0.25)**(2./3.)-0.2241*SIN(WRK5)/(K*(AI+0.25)**(2 1./3.)) 5 CONTINUE 3 CONTINUE DO 6 I=1,15 6 TAU(I)=CMPLX(TAUR(I),TAUI(I)) DO 7 I=1,15 7 E(I)=(752./D)*SQRT(CHI)*CEXP(J*TAU(I)*CHI)/(2.*TAU(I)-1./DELTA**2) E2=CABS(E(1)+E(2)+E(3)+E(4)+E(5)+E(6)+E(7)+E(8)+E(9)+E(10)+E(11)+E 1(12)+E(13)+E(14)+E(15)) IF(E2.LT.1.E-06) GO TO 8 DB=20.*ALOG10(1000.*E2) 8 CONTINUE reTURN END