standev2.f

C
C        PROGRAM STANDEV2.F 
C
C     THIS ROUTINE PLOTS THE ROOT MEAN SQUARE OF THE DIFFERENCE BETWEEN
C     THE LOG OF THE OBSERVED RECORD SPECTRA (D.DAT) 
C     FROM PROGRAM QUARTET.F AND THE LOG OF THE MODEL ESTIMATED RECORD SPECTRA 
C     FROM PROGRAM RUNSPARS1.F.
C
C     ALSO WRITES OUT A FILE OF THE NUMBER OF SOURCES FOR EACH SITE,
C     NUMPERSITE.DAT
C
      DIMENSION D(1000),DFIT(200000),DF(1000),FREQ(1000)
      DIMENSION H(100),XMAG(100),STALAT(100),STALON(100),ANGACC(100),
     +ANGVEL(100),IASF(13),WINDOW(4),NUM(100),
     +DATA(100,100,100),SD(100,100),SYN(100,100,100)
      DIMENSION SORLAT(100),SORLON(100),CHANALT(3,100),CHANAZI(3,100)
      CHARACTER*5 BLANK
      CHARACTER*7 BLANK2
      CHARACTER*40 STRING,STALOC,AFLOC
      CHARACTER*14 AFILE(100),NAME1
      CHARACTER*40 TITLE,NAME
      CHARACTER*5 STA(100),NAME2
      DATA IASF /13*1/
      CALL SETUSV('MU',8)
      CALL OPNGKS
      OPEN(UNIT=15,FILE='NUMPERSITE.DAT',STATUS='NEW')
      OPEN(UNIT=11,FILE='STANDEV2.DAT',STATUS='NEW')
      WRITE(6,660)
 660  FORMAT(1X,'WHAT LINE WEIGHT THICKNESS DO YOU WANT',/,
     +1X,'THE DEFAULT VALUE IS 1.')
      READ(5,*) WTLINE
      CALL GSLWSC(WTLINE)
      CALL GSASF(IASF)
      CALL GSCR(1,0,0.,0.,1.)
      CALL GSCR(1,1,1.,1.,1.)
      CALL FRAME
      CALL SET(0.,1.,0.,1.,0.,1.,0.,1.,1)
      ICOUNT=1
      WRITE(6,600)
 600  FORMAT(1X,'HOW MANY DATA RECORDS ARE THERE')
      READ(5,*) ICREC
      OPEN(UNIT=7,FILE='D.DAT',STATUS='OLD',ACCESS='SEQUENTIAL',
     +FORM='UNFORMATTED')
      OPEN(UNIT=9,FILE='DFIT.DAT',STATUS='OLD',ACCESS='SEQUENTIAL',
     +FORM='UNFORMATTED')
      READ(9) NRA
      READ(9) (DFIT(I),I=1,NRA)
      CLOSE(9)
C
C     TAKE INVERSE LOGARITHM
C
C      DO 21 I=1,NRA
C  21  DFIT(I)=10.**DFIT(I)
C
C     ZERO OUT ARRAYS.
C
      DO 24 I=1,100
      DO 24 J=1,100
  24  SD(J,I)=0.
      DO 25 I=1,100
  25  NUM(I)=0
C
      WRITE(6,621)
 621  FORMAT(1X,'SPECIFY BACKGROUND TYPE FOR PLOTS:',/,
     +1X,'1 = PERIMETER BACKGROUND',/,
     +1X,'2 = GRID BACKGROUND',/,
     +1X,'3 = JUST AXES',/,
     +1X,'4 = NO BACKGROUND')
      READ(5,*) IBAK
      WRITE(6,622)
 622  FORMAT(1X,'SPECIFY LINEAR (=0) OR LOGARITHMIC (=1) X, HORIZ.',/,
     +1X,'SCALE')
      READ(5,*) IXTYPE
      WRITE(6,623)
 623  FORMAT(1X,'SPECIFY LINEAR (=0) OR LOGARITHMIC (=1) Y, VERT.',/,
     +1X,'SCALE')
      READ(5,*) IYTYPE
C
C     FIRST FIND OUT HOW MANY DIFFERENT SOURCES THERE ARE.
C
      WRITE(6,707)
 707  FORMAT(1X,'ENTER THE NAME OF THE FILE IN MASTER.DAT FORMAT',/,
     +1X,'WITH SOURCE LOCATIONS')
      READ(5,601) AFLOC
      OPEN(UNIT=10,FILE=AFLOC,STATUS='OLD')
      IC=1
  4   READ(10,602,END=5) XLAT,XLONG,H(IC),XMAG(IC),AFILE(IC) 
 602  FORMAT(F7.4,1X,F9.4,2F5.2,1X,A)
      SORLAT(IC)=XLAT
      SORLON(IC)=XLONG
      IC=IC+1
      GO TO 4
  5   IC=IC-1
      CLOSE(10)
      NSOR=IC
      WRITE(6,666) (SORLAT(I),SORLON(I),I=1,NSOR)
 666  FORMAT(1X,2F12.3)     
C
C     NOW READ IN SITE INFORMATION.
C
      WRITE(6,702)
 702  FORMAT('ENTER THE NAME OF THE FILE WITH STATION LOCATIONS')
      READ(5,601) STALOC
      OPEN(UNIT=10,FILE=STALOC,STATUS='OLD')
      I=1
  61  READ(10,601,END=60) STRING
      IF(STRING(1:5) .EQ. 'STNID') THEN 
      STA(I)=STRING(7:11)
      READ(10,*) BLANK,STALAT(I),STALON(I),HH
      DO 37 IPC=1,3
      READ(10,*,END=60) BLANK2,ICOM,CHANALT(ICOM,I)
      READ(10,*,END=60) BLANK2,ICOM,CHANAZI(ICOM,I)
  37  CONTINUE
      I=I+1
      END IF
      GO TO 61
  60  CLOSE(10)
 601  FORMAT(A)
      NSIT=I-1
C
C
      DO 100 IK=1,ICREC
      READ(7) ILOOP,JLOOP,IFN,FMIN,FSTEP
      K1=IFN*(IK-1)+1
      K2=IFN*IK
      READ(7) (D(K),K=1,IFN)
C
C     TAKE INVERSE LOGARITHM
C
C      DO 22 I=1,IFN
C  22  D(I)=10.**D(I)
C
      NAME1=AFILE(ILOOP)
      NAME2=STA(JLOOP)
      TITLE(1:33)='EVENT: '//NAME1(1:14)//' SITE: '//NAME2(1:5)
      IF(IK .EQ. 1) THEN
      DO 20 I=1,IFN
  20  FREQ(I)=FMIN+FLOAT(I-1)*FSTEP
      END IF
      KQ=K1
      DO 10 I=1,IFN
      DF(I)=DFIT(KQ)
  10  KQ=KQ+1
      NUM(JLOOP)=NUM(JLOOP)+1
      DO 23 I=1,IFN
      SYN(I,JLOOP,NUM(JLOOP))=DF(I)
  23  DATA(I,JLOOP,NUM(JLOOP))=D(I)
 100  CONTINUE
      WRITE(6,28)
      WRITE(6,27) (STA(I),NUM(I),I=1,NSIT)
      WRITE(15,27) (STA(I),NUM(I),I=1,NSIT)
      CLOSE(15)
  27  FORMAT(1X,A,5X,I5)
  28  FORMAT(1X,'NUMBER OF SOURCES FOR EACH SITE')
      DO 31 I=1,NSIT
      DO 31 IR=1,NUM(I)
      DO 31 J=1,IFN
  31  SD(J,I)=SD(J,I)+(DATA(J,I,IR)-SYN(J,I,IR))**2
      DO 33 I=1,NSIT
C*****This next term is for RMS not standard deviation****
      NV=NUM(I)
      IF(NV .LT. 1) THEN
      DO 39 J=1,IFN
  39  SD(J,I)=0.
      END IF
      IF(NV .GT. 0) THEN
      DO 38 J=1,IFN
  38  SD(J,I)=10.**SQRT(SD(J,I)/FLOAT(NV))
      END IF
  33  CONTINUE
      DO 200 IMQ=1,NSIT
      DO 29 I=1,IFN
  29  DF(I)=SD(I,IMQ)
      WRITE(11,603) STA(IMQ)
      WRITE(11,604) (DF(I),I=1,IFN)
 604  FORMAT(1X,8F12.5)
 603  FORMAT(A)
      NAME(1:28)='STANDARD DEV. AT SITE: '//STA(IMQ)(1:5)
      CALL BOX4(FREQ,DF,IFN,WINDOW,IBAK,IXTYPE,IYTYPE,NAME,
     +ICOUNT)
 200  CONTINUE
      CLOSE(7)
      CLOSE(11)
      CALL CLSGKS
      STOP
      END
      SUBROUTINE DRAW(X,Y,NT,WINDOW,IBAK,IXTYPE,IYTYPE,NAME)
C
C     X......ARRAY OF X (HORIZONTAL) COORDINATE VALUES
C     Y......ARRAY OF Y (VERTICAL) COORDINATE VALUES
C     NT.....NUMBER OF POINTS IN X AND Y
C     WINDOW.FOUR ELEMENT ARRAY SPECIFYING THE GRAPH WINDOW
C            AS FRACTIONS OF THE SCREEN DIMENSIONS.
C            WINDOW(1)= LEFT LIMIT
C            WINDOW(2)= RIGHT LIMIT
C            WINDOW(3)= BOTTOM LIMIT
C            WINDOW(4)= TOP LIMIT.
C     IBAK...BACKGROUND PARAMETER.
C            IBAK= 1, PERIMETER BACKGROUND
C            IBAK= 2, GRID BACKGROUND
C            IBAK= 3, JUST AXES
C            IBAK= 4, NO BACKGROUND
C     IXTYPE.SPECIFIES LINEAR (=0) OR LOGARITHMIC (=1) X SCALE.
C     IYTYPE.SPECIFIES LINEAR (=0) OR LOGARITHMIC (=1) Y SCALE.
C     NAME...TITLE FOR PLOT.
C
      DIMENSION X(1),Y(1),WINDOW(4)
      CHARACTER*12 LNLG(2)
      CHARACTER*40 NAME
      IF(IXTYPE .EQ. 0) LNLG(1)='LINEAR$'
      IF(IXTYPE .EQ. 1) LNLG(1)='LOGARITHMIC$'
      IF(IYTYPE .EQ. 0) LNLG(2)='LINEAR$'
      IF(IYTYPE .EQ. 1) LNLG(2)='LOGARITHMIC$'
C      CALL AGSETF('Y/MINIMUM.',0.)
C      CALL AGSETF('Y/MAXIMUM.',2.)
C      CALL AGSETF('Y/SMALLEST.',0.)
C      CALL AGSETF('Y/LARGEST.',2.)
      CALL AGSETI('FRAME.',2)
      CALL AGSETP('GRAPH WINDOW.',WINDOW(1),4)
      CALL AGSETI('BACKGROUND TYPE.',IBAK)
      CALL AGSETI('X/LOGARITHMIC.',IXTYPE)
      CALL AGSETI('Y/LOGARITHMIC.',IYTYPE)
      CALL AGSETC('LABEL/NAME.','B')
      CALL AGSETI('LINE/NUMBER.',-100)
      CALL AGSETC('LINE/TEXT.',LNLG(1))
      CALL AGSETC('LABEL/NAME.','L')
      CALL AGSETI('LINE/NUMBER.',100)
      CALL AGSETC('LINE/TEXT.',LNLG(2))
C      DO 62 I=1,NT
C      IF(Y(I) .LT. 0.) Y(I)=0.
C  62  IF(Y(I) .GT. 2.) Y(I)=2.
      CALL EZXY(X,Y,NT,NAME)
      CALL SET(0.,1.,0.,1.,0.,1.,0.,1.,1)
      RETURN
      END
      SUBROUTINE BOX4(X,Y,NT,WINDOW,IBAK,IXTYPE,IYTYPE,NAME,
     +ICOUNT)
      DIMENSION X(1),Y(1),WINDOW(4)
      CHARACTER*40 NAME
      IF(ICOUNT .GT. 4) THEN
      CALL FRAME
      ICOUNT=1
      END IF
      IF(ICOUNT .EQ. 1) THEN
      WINDOW(1)=0.
      WINDOW(2)=0.5
      WINDOW(3)=0.5
      WINDOW(4)=1.0
      END IF
      IF(ICOUNT .EQ. 2) THEN
      WINDOW(1)=0.5
      WINDOW(2)=1.0
      WINDOW(3)=0.5
      WINDOW(4)=1.0
      END IF
      IF(ICOUNT .EQ. 3) THEN
      WINDOW(1)=0.
      WINDOW(2)=0.5
      WINDOW(3)=0.
      WINDOW(4)=0.5
      END IF
      IF(ICOUNT .EQ. 4) THEN
      WINDOW(1)=0.5
      WINDOW(2)=1.0
      WINDOW(3)=0.
      WINDOW(4)=0.5
      END IF
      CALL DRAW(X,Y,NT,WINDOW,IBAK,IXTYPE,IYTYPE,NAME)
      ICOUNT=ICOUNT+1
      RETURN
      END