www.pudn.com > Noise.rar > NOISBW.FOR
PROGRAM NOISBW
C.....A FORTRAN V ROUTINE FOR PC USING USER INPUT OF SEASON, LOCATION,
C.....FREQUENCY, BANDWIDTH AND LOCAL TIME AND 8 EXTERNAL FILES OF COEFFICIENTS
C.....TO GIVE OUTPUT OF ATMOSPHERIC NOISE WITH STATISTICAL VALUES FOR DL, DU,
C.....SL, SM, SU, VD, AND SVD. THIS IS A MODIFICATION (1/89) OF THE
C.....PROGRAM NOISB.
COMMON /ANOIS/ ATNU,ATNY,CC,TM,RCNSE,DU,DL,SIGM,SIGU,SIGL,KJ,JK
COMMON /TON/ ATMO, GNOS, ZCNSE, XADJN, XNOISE, ZNOISE
COMMON /TWO/ DUD(5,12,5),FAM(14,12),FAKP(29,16,6),FAKABP(2,6)
COMMON /THREE/ VDARRAY(5,12,2)
COMMON /NSTAT/ DLA,DUA,SLA,SMA,SUA,VD200,SVD200
DIMENSION FREQA(60),TBHR(6),FREQL(11)
DIMENSION AARR(4,60,6),VDARR(4,60,6),SVDARR(4,60,6),
A DLARR(4,60,6),DUARR(4,60,6),SLARR(4,60,6),
B SMARR(4,60,6),SUARR(4,60,6)
CHARACTER*6 SEASON(4),SEAFIN(4),VDFIN(4)
CHARACTER*9 TIMEBLK(6)
CHARACTER*40 LOCNAME
DATA SEASON /'WINTER','SPRING','SUMMER','AUTUMN'/
DATA SEAFIN /'ISCOF1','ISCOF2','ISCOF3','ISCOF4'/
DATA VDFIN /'VDCOF1','VDCOF2','VDCOF3','VDCOF4'/
DATA TIMEBLK /'0000-0400','0400-0800','0800-1200','1200-1600',
A '1600-2000','2000-2400'/
DATA TBHR /2.0,6.0,10.0,14.0,18.0,22.0/
DATA FREQL /.01,.02,.05,.1,.2,.5,1.,2.,5.,10.,20./
C.....
C.....BEGIN INPUT WITH SEASON
C.....
100 WRITE (*,1050)
1050 FORMAT (/1X,' INPUT SEASON, 1=WINTER, 2=SPRING, 3=SUMMER'/
A ' 4=AUTUMN, 5=ALL SEASONS, ANYTHING ELSE=END OF RUN')
READ (*,*) ISEASON
ISEAIN=ISEASON
IF (ISEASON .LT. 1 .OR. ISEASON .GT. 5) GO TO 500
ISE1=ISEASON
ISE2=ISEASON
IF (ISEASON .EQ. 5) ISE1=1
IF (ISEASON .EQ. 5) ISE2=4
DO 275 ISE = ISE1,ISE2
ISEASON = ISE
IF (ISEAIN .GT. 4 .AND. ISE .GT. 1) GO TO 105
C.....
C.....INPUT LOCATION
C.....
WRITE (*,1051)
1051 FORMAT (/1X,' INPUT LOCATION LATITUDE (- IF S), LONGITUDE (- IF'
A ' WEST)'/' AND LOCATION NAME IN SINGLE QUOTES')
READ (*,*) RLATD,ALONGD,LOCNAME
RLONGD=ALONGD
IF (ALONGD .LT. 0.0) RLONGD=360.+ALONGD
105 CONTINUE
KODESEA=ISEASON
IF (RLATD .GE. 0.0) GO TO 110
IF (ISEASON .EQ. 1) KODESEA = 3
IF (ISEASON .EQ. 2) KODESEA = 4
IF (ISEASON .EQ. 3) KODESEA = 1
IF (ISEASON .EQ. 4) KODESEA = 2
110 CONTINUE
OPEN (3,FILE=SEAFIN(KODESEA),STATUS='UNKNOWN',ACCESS=
A 'SEQUENTIAL', FORM='BINARY',MODE='READ')
READ (3,END=300,ERR=400) DUD,FAM,FAKP,FAKABP
CLOSE (3)
OPEN (5,FILE=VDFIN(KODESEA),STATUS='UNKNOWN',ACCESS=
A 'SEQUENTIAL', FORM='BINARY',MODE='READ')
READ (5,END=300,ERR=400) VDARRAY
CLOSE (5)
LUFO=4
OPEN (4,FILE='NOISBW.OUT',STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
A FORM='FORMATTED')
IF (ISE1 .NE. ISE2 .AND. ISE .NE. 1) GO TO 135
C.....
C.....INPUT BANDWIDTH
C.....
WRITE (*,1052)
1052 FORMAT (/1X,' INPUT BANDWIDTH'/
C ' AND OUTPUT CODE 1 FOR DBW OR 2 FOR FA')
READ (*,*) BW,IODBWDB
BWR=BW/200.
C.....
C.....INPUT FMHZ
C.....
WRITE (*,1053)
1053 FORMAT (/1X,' INPUT SPECIFIC FMHZ DESIRED (.01-30.)'/
A ' OR 0.0 FOR COMPUTER ARRAY LOOP OF FMHZ'/
B ' OR -1.0 FOR USER ARRAY LOOP OF FMHZ'/
C ' OR -2.0 FOR LOGARITHMIC ARRAY LOOP OF FMHZ.'/
D ' DO NOT ASK FOR MORE THAN 60 VALUES OF FMHZ.')
READ (*,*) FREQ
NF=1
FREQA(NF)=FREQ
IF (FREQ-0.0) 125,122,130
122 WRITE (*,1054)
1054 FORMAT (/1X,' INPUT BEGINNING, ENDING, AND INCREMENT FOR FREQ')
READ (*,*) FBEG,FEND,FINC
FREQA(NF)=FBEG
IF (FBEG .LT. .01 .OR. FEND .GT. 30.) GO TO 295
DO 124 KF=2,60
NF=NF+1
FREQA(NF)=FREQA(NF-1)+FINC
IF ((FREQA(NF)+.005)-FEND) 124,129,129
124 CONTINUE
GO TO 295
125 IF (FREQ .NE. -1.0) GO TO 127
DO 126 KF=1,61
WRITE (*,1057)
1057 FORMAT (/1X,' INPUT FREQUENCY (0.0 TO END)')
READ (*,*) FREQA(NF)
IF (FREQA(NF) .EQ. 0.0) GO TO 129
NF=NF+1
126 CONTINUE
GO TO 295
127 IF (FREQ .NE. -2.0) GO TO 295
NF=11
DO 128 KF=1,11
FREQA(KF)=FREQL(KF)
128 CONTINUE
GO TO 130
129 IF ((FREQ .EQ. -1.) .OR. (FREQ .EQ. 0. .AND. FREQA(NF) .GT.
A FEND)) NF=NF-1
C.....
C.....INPUT TIME
C.....
130 WRITE (*,1055)
1055 FORMAT (/1X,' INPUT SPECIFIC LOCAL MEAN TIME (0-23)'/
A ' OR 25 FOR TIME BLOCKS')
READ (*,*) RLMT
IF (RLMT .LT. 0. .OR. RLMT .EQ. 24. .OR. RLMT .GT. 25.)
A GO TO 295
135 CONTINUE
IF (RLMT-25.) 140,200,295
C.....
C.....USE THIS BRANCH TO DO SPECIFIC TIME
C.....
140 CALL ANOIS1(RLMT,RLATD,RLONGD)
WRITE (LUFO,1002) RLATD,ALONGD,LOCNAME
1002 FORMAT(////2X,'LAT = ',F6.2,', LONG = ',F7.2,', ',A40)
WRITE (LUFO,1014) SEASON(ISEASON),RLMT,BW
1014 FORMAT (2X,A6,', LMT = ',F4.1,', BANDWIDTH = ',F7.0/)
IF (IODBWDB .EQ. 1) WRITE (LUFO,1005)
IF (IODBWDB .EQ. 2) WRITE (LUFO,1013)
1005 FORMAT (13X,'---MEDIAN ATMOSPHERIC NOISE OR DBW(1HZ)--')
1013 FORMAT (13X,'--MEDIAN ATMOSPHERIC NOISE, FA (DB>KTO)--')
WRITE (LUFO,1008)
1008 FORMAT (7X,'FMHZ',3X,'NOISE',
C 3X,'DL',3X,'DU',3X,'SL',3X,'SM',3X,'SU'3X,'VD',2X,'SVD'/)
DO 150 IFREQ=1,NF
FREQ=FREQA(IFREQ)
CALL MGENOIS(FREQ,RLATD)
IF (IODBWDB .EQ. 1) GO TO 148
ATMO=ATMO+204.
148 CONTINUE
VD=VDC(VD200,BWR)
SVD=SVDC(VD200,SVD200,BWR)
WRITE (LUFO,1003) FREQ,ATMO,
A DLA,DUA,SLA,SMA,SUA,VD,SVD
1003 FORMAT(4X,F7.3,F8.1,7F5.1)
150 CONTINUE
WRITE (LUFO, 1018)
1018 FORMAT (////)
GO TO 275
C.....
C.....USE THIS BRANCH TO DO TIME BLOCKS
C.....
200 CONTINUE
DO 270 IFREQ=1,NF
FREQ=FREQA(IFREQ)
DO 265 ITB=1,6
TIME=TBHR(ITB)
CALL ANOIS1(TIME,RLATD,RLONGD)
CALL MGENOIS(FREQ,RLATD)
IF (IODBWDB .EQ. 1) GO TO 258
ATMO=ATMO+204.
258 CONTINUE
AARR(ISE,IFREQ,ITB)=ATMO
DLARR(ISE,IFREQ,ITB)=DLA
DUARR(ISE,IFREQ,ITB)=DUA
SLARR(ISE,IFREQ,ITB)=SLA
SMARR(ISE,IFREQ,ITB)=SMA
SUARR(ISE,IFREQ,ITB)=SUA
VDARR(ISE,IFREQ,ITB)=VDC(VD200,BWR)
SVDARR(ISE,IFREQ,ITB)=SVDC(VD200,SVD200,BWR)
265 CONTINUE
270 CONTINUE
275 CONTINUE
C.....
C.....OUTPUT FOR TIME BLOCK COMPUTATIONS; ELSE GO FOR NEW INPUT
C.....
IF (RLMT-25.) 100,280,295
280 CONTINUE
DO 286 IFREQ=1,NF
DO 284 ISE=ISE1,ISE2
WRITE (LUFO,1002) RLATD,ALONGD,LOCNAME
WRITE (LUFO,1004) SEASON(ISE),FREQA(IFREQ),BW
1004 FORMAT (2X,A6,', FMHZ = ',F6.3,', BANDWIDTH = ',F7.0/)
IF (IODBWDB .EQ. 1) WRITE (LUFO,1005)
IF (IODBWDB .EQ. 2) WRITE (LUFO,1013)
WRITE (LUFO,1006)
1006 FORMAT (1X,'TIME BLOCK',3X,'NOISE',
C 3X,'DL',3X,'DU',3X,'SL',3X,'SM',3X,'SU',3X,'VD',2X,'SVD'/)
DO 282 ITB = 1,6
WRITE (LUFO,1007) TIMEBLK(ITB),AARR(ISE,IFREQ,ITB),
C DLARR(ISE,IFREQ,ITB),DUARR(ISE,IFREQ,ITB),SLARR(ISE,IFREQ,ITB),
C SMARR(ISE,IFREQ,ITB),SUARR(ISE,IFREQ,ITB),VDARR(ISE,IFREQ,ITB),
C SVDARR(ISE,IFREQ,ITB)
1007 FORMAT (2X,A9,F8.1,7F5.1)
282 CONTINUE
284 CONTINUE
WRITE (LUFO,1018)
286 CONTINUE
GO TO 100
295 WRITE (*,2010)
GO TO 100
300 WRITE (*,2020)
CLOSE (4)
STOP
400 WRITE (*,2030)
CLOSE (4)
STOP
500 CONTINUE
CLOSE (4)
WRITE (*,2000)
STOP
2000 FORMAT (//' NORMAL PROGRAM TERMINATION, OUTPUT ON NOISBW.OUT')
2010 FORMAT (//' ERROR ON INPUT, TRY AGAIN')
2020 FORMAT (//' END OF FILE ON DATA BASE')
2030 FORMAT (//' ERROR ON DATA BASE READ')
END
SUBROUTINE ANOIS1(RLMT,RLATD,RLONGD)
CR....A ROUTINE THAT USES RLMT TO DETERMINE THE TIMEBLOCK (KJ)
CR....AND ADJACENT TIME BLOCK (JK) (THIS IS THE PRIOR TIMEBLOCK
CR....FOR THE FIRST 2 HOURS OF KJ, THE SAME, IE JK=KJ, FOR THE 3RD
CR....HOUR OF KJ AND THE NEXT TIME BLOCK FOR THE LAST HOUR OF KJ)
CR....AND THEN CALLS NOISY TO FIGURE THE ATMOSPHERIC NOISE (ATNU
CR....OR ATNY) FOR EACH OF THESE TIME BLOCKS.
C.....
C.....THIS ROUTINE DETERMINES THE 1 MHZ ATMOSPHERIC NOISE
C.....
C.....FOURIER SERIES IN LATITUDE AND LONGITUDE FOR TWO DISCRETE
C.....LOCAL TIME BLOCKS
C.....
COMMON /ANOIS/ ATNU,ATNY,CC,TM,RCNSE,DU,DL,SIGM,SIGU,SIGL,KJ,JK
C.....LMT AT RCVR SITE
100 CC = RLMT
KJ= 6
IF(CC-22.) 105,110,110
105 KJ = CC/4. +1.
110 TM = 4*KJ-2
IF(CC-TM) 115,120,125
115 JK = KJ -1
GO TO 130
120 JK = KJ
GO TO 130
125 JK = KJ+1
130 IF(JK) 135,135,140
135 JK =6
GO TO 150
140 IF(JK-6) 150,150,145
145 JK = 1
C.....EAST LONGITUDE (IN DEGREES)
150 CEG= RLONGD
165 XLA =RLATD
C.....LATITUDE (IN DEGREES) "+" IS NORTH
CALL NOISY(KJ,XLA,CEG,ATNU)
CALL NOISY(JK,XLA,CEG,ATNY)
RETURN
END
SUBROUTINE NOISY (KJ, XLA, CEG, ANOS)
CR....A ROUTINE TO USE THE TIMEBLOCK (KJ), THE LAT (XLA), THE LONG
CR....(CEG), AND THE COEFFICIENTS (FAKP AND FAKAB) TO
CR....DETERMINE THE ATMOSPHERIC NOISE (ANOS).
CR....THIS ROUTINE USES MAPS TO GET THE 1 MHZ FAM VALUE.
C.....NOISY IS A GENERAL PURPOSE ROUTINE USED TO EVALUATE A FOURIER
C.....SERIES IN TWO VARIABLES.
C.....KJ --- NUMBER OF FOURIER COEFFICIENT ARRAY TO BE USED
C.....XLA --- GEOGRAPHIC LATITUDE, DEGREES,
C.....CEG --- GEOGRAPHIC EAST LONGITUDE, DEGREES
C.....ANOS --- NOISE VALUE, MEDIAN POWER DB ABOVE KTB
C.....FAKABP --- NORMALIZING FACTORS FOR FOURIER SERIES
C.....KJ = 1 TO 6 IS ATMOSPHERIC NOISE
C.....
C.....* NOTE - XLA, CEG, ANOS, FAKABP ARE NOT ALWAYS AS PREVIOUSLY
C..... DEFINED
C.....FOURIER VARIABLES AND ATMOSPHERIC RADIO NOISE
C.....
COMMON /TWO/ DUD(5,12,5),FAM(14,12),FAKP(29,16,6),FAKABP(2,6)
DIMENSION SX (15), SY(29), ZZ (29)
IF (KJ - 6) 105,105,100
100 KJ=6
C.....LIMITS OF FOURIER SERIES
105 LM = 29
LN = 15
C.....HALF ANGLE (IN RADIANS)
110 Q = .0087266466 * CEG
C.....LONGITUDE SINES
DO 115 K = 1, 15
115 SX(K)=SIN(Q*K)
118 CONTINUE
C.....LONGITUDE SERIES
DO 125 J = 1, LM
R = 0.
DO 120 K = 1, LN
120 R = R + SX (K) * FAKP (J, K, KJ)
125 ZZ (J) = R + FAKP (J, 16, KJ)
C.....ANGLE PLUS 90 DEGREES (IN RADIANS)
Q = .01745329252 * (XLA + 90.)
C.....LATITUDE SERIES
DO 140 J=1,29
140 SY(J)=SIN(Q*J)
R = 0.
DO 130 K = 1, LM
130 R = R + SY (K) * ZZ (K)
C.....FINAL FOURIER SERIES EVALUATION (NOTE LINEAR NORMALIZATION)
135 ANOS = R + FAKABP(1,KJ)+FAKABP(2,KJ)* Q
RETURN
END
SUBROUTINE MGENOIS(FREQ,RLAT)
C.....
CR....THIS ROUTINE IS A MODIFIED VERSION OF GENOIS WHERE JUST THE ATMOSPHERIC
CR....NOISE STATISTICS ARE CALCULATED AND VD AND SVD ARE ADDED.
C.....
COMMON /ANOIS/ ATNU,ATNY,CC,TM,RCNSE,DU,DL,SIGM,SIGU,SIGL,KJ,JK
COMMON /TON/ ATMO, GNOS, ZCNSE, XADJN, XNOISE, ZNOISE
COMMON /NSTAT/ DLA,DUA,SLA,SMA,SUA,VD200,SVD200
DUME = AMIN1(FREQ,55.)
CALL MGENFAM(RLAT,KJ,DUME,ATNU,ATNZ,DU,DL,SIGM,SIGU,SIGL,VD1,SVD1)
CALL MGENFAM(RLAT,JK,DUME,ATNY,ATNX,DX,DQ,SIGZ,SIGX,SIGSQ,VD2,
C SVD2)
C.....BEGIN OF INTERPOLATION ON LOCAL TIME
SLOP = ABS(CC-TM)/4.
ATNOS = ATNZ + (ATNX - ATNZ) * SLOP
ATMO=ATNOS-204.
DUA= DU +(DX-DU)*SLOP
DLA= DL +(DQ-DL)*SLOP
SMA= SIGM+ (SIGZ-SIGM)*SLOP
SUA= SIGU +(SIGX-SIGU)*SLOP
SLA= SIGL+(SIGSQ-SIGL)* SLOP
VD200=VD1+(VD2-VD1)*SLOP
SVD200=SVD1+(SVD2-SVD1)*SLOP
RETURN
END
SUBROUTINE MGENFAM(Y2,IBLK,FREQ,Z,FA,DU,DL,DMS,DUS,DLS,VD,SVD)
c**********************************************************************
c Re-written 3.June.93 by Greg Hand because previous version was
c really incorrect. It made an attempt to limit Sigma Fam (DMS)
c to a 10 MHz frequency, but the indicies I and J became
c confused, and the result was not correct. This current
c version should limit DMS to 10 MHz and the others to 20 MHz
c because the curves end at 20 MHz. Unfortunately, this error
c has probably existed since time began, and it may take a
c while for this corrected version to propagate into all version
c that exist. The magnitude of the error that would have been
c caused is not known, but it is believed to be small.
c**********************************************************************
C.....
CR....THIS IS A MODIFIED GENFAM ROUTINE; AN ADDITIONAL SECTION DEFINES VD, SVD.
C.....GENFAM CALCULATES THE FREQUENCY DEPENDENCE OF THE ATMOSPHERIC
C.....NOISE AND GETS DECILES AND PREDICTION ERRORS FROM TABLES
C.....
COMMON /TWO/ DUD(5,12,5),FAM(14,12),FAKP(29,16,6),FAKABP(2,6)
COMMON /THREE/ VDARRAY(5,12,2)
DIMENSION V(5)
IBK=IBLK
C.....CHECK IF LATITUDE IS NORTH OR SOUTH
IF (Y2.lt.0.) ibk=ibk+6
U1 = - .75
X = ALOG10(FREQ)
U = (8. * 2. * * X - 11.) / 4.
KOP = 1
110 PZ = U1 * FAM (1, IBK) + FAM (2, IBK)
PX = U1 * FAM (8, IBK) + FAM (9, IBK)
DO 115 I = 3, 7
PZ = U1 * PZ + FAM (I, IBK)
115 PX = U1 * PX + FAM (I + 7, IBK)
IF(KOP.eq.1) then
CZ = Z * PZ + PX
CZ = Z + Z - CZ
U1 = U
KOP = 2
GO TO 110
end if
FA = CZ * PZ + PX
c Limit frequency to 20 MHz for DUA, DLA, DUS, DLS
c because curves in REP 322 only go to 20 MHz
if(FREQ.gt.20.) X=ALOG10(20.)
DO 145 I = 1, 5
c Limit frequency to 10 MHz for DMS (Sigma Fam)
c because curves in REP 322 only go to 10 MHz
if(I.eq.5 .and. FREQ.gt.10.) X=1.
Y = DUD(1,IBK,I)
DO 140 J = 2,5
140 Y = Y*X + DUD(J,IBK,I)
145 V(I) = Y
DU = V (1)
DL = V (2)
DUS = V (3)
DLS = V (4)
DMS = V (5)
X = ALOG10(FREQ)
if(FREQ.gt.20.) X=ALOG10(20.)
DO 165 I = 1, 2
Y = VDARRAY (1, IBK, I)
do 160 j=2,5
160 Y = Y * X + VDARRAY (J, IBK, I)
165 V (I) = Y
VD = V (1)
SVD = V (2)
RETURN
END
FUNCTION VDC(VD200,BWR)
C OBTAINS THE NOISE PARAMETER -VD-, FOR THE SPECIFIED BANDWIDTH
C FROM THE CCIR REPORT 322 (OR OTHER) 200HZ BANDWIDTH -VD- (VD200),
C -BWR- IS THE BANDWIDTH RATIO (REQUIRED BANDWIDTH/200 HZ BANDWIDTH).
IF (VD200 .LE. 1.049) GO TO 50
VDO=VD200+(0.4679+0.2111*VD200)*ALOG10(BWR)
IF (VDO .LE. 1.049) GO TO 50
VDC=VDO
RETURN
50 VDC=1.049
RETURN
END
FUNCTION SVDC(VD200,SVD200,BWR)
V01=VDC((VD200+SVD200),BWR)
V02=VDC((VD200-SVD200),BWR)
SVDC=(V01-V02)/2.
RETURN
END