C spectres a partir d'un fichier ah (un seul fichier par evenement)
c    14 mars 1990 B.R.
C  CALCUL DES SPECTRES DE TRAINS D'ONDES DE RAYLEIGH
C       HYPO=0. POUR STATIONS GEOSCOPE AVEC POLARITE NORMALE
C       HYPO=PI  SINON
	include 'ahhead_common.bis'
        character*8 idis,ipl
       character*80 inahfile,filout
       character*6 stcod,chan,ist,ichan
       character*1 tr,tq
      DIMENSION TT(500),TB(500),ichan(5) 
      INTEGER*2 INM(90,4),IND(90,5)
      INTEGER*2 ISTA(5)
      INTEGER*2 IFIL(12),fipl(4) 
      COMMON/XRR/Y(20000),X(200000),YY(30000),T(30000),XX(500)
      COMMON/HUGE/XW(200000),U1(4),U2(4),U10(4),U20(4)
      COMMON/HEADIR/A0,POLES(30,4),ZEROS(30,4),coen(60,5),coed(60,5),
     *dtp(4),dtc(5),nsp,npp(4),nzz(4),nsc,ncn(5),ncd(5),kuni,tr(4),tq(5)
	common/polzer/fmul,zer(10),pol(20),nzi,npi,fnorm,nfilts,
     *nftyp(6),nlen(6),dti(6),ci(500)
c      INTEGER*2 ANON,LH,LM,LFH,LFM,L,CN
      INTEGER*2 ANON,CN
c      INTEGER*2 INAM(5),ITITRE(10),IH,IM,LA,LMO,LJ
      INTEGER*2 INAM(5),ITITRE(10) 
      COMPLEX C(12000),pole(30),zero(20),pole1(30),zero1(30),r,poles,zeros
	complex pol,zer,evres
      EQUIVALENCE(IFIL,IDIS)
      EQUIVALENCE(ISTA,PAREIL)
      equivalence(fipl,ipl)
      DATA ANON/'NO'/
      DATA CN/'N'/,NKO/' '/
      DATA INAM/'1','2','3','4','5'/
      DATA JKO,JKI/'  ','  '/
      DATA IFIL(1)/'DI'/,IFIL(2)/'ST'/,fipl/'pl','sp','da','ta'/
       print*,'creat distance filename'
      READ(5,99)(IFIL(I),I=3,12)
      write(6,99)(IFIL(I),I=3,12)
99     FORMAT(10A2)
      print*,'open idis'
      OPEN(12,FILE=IDIS)
      read(5,"(a80)")filout 
      write(6,"(a80)")filout
      print*,'open filout'
       open(9,file=filout)
       read(5,"(a80)")inahfile
       write(6,"(a80)")inahfile
       luin=15
      open(11,file=ipl)
c     PRINT*,'DATE OF EVENT'
c      READ(5,120)LA,LMO,LJ
c      PRINT*,'ORIGIN TIME OF EVENT'
c      READ(5,122)IH1,IM,SEC
122      FORMAT(I2,I3,F5.1)
120    FORMAT(I2,1X,I2,1X,I2)
C  TT1.GT.TT2 LIMITES DES PERIODES CONSERVEES
C  EN GENERAL 400 ET 100 SEC.
      PRINT*,'PERIODS OF SPECTRUM'
      READ(5,*)TT1,TT2
       print*,TT1,TT2
      READ(5,*)(U10(I),U20(I),I=1,4)
      print*,(U10(I),U20(I),I=1,4)
      PI=3.1415926535
      TPI=2.*PI
      PRINT*,'WRITE TIME SERIES YES=1,NO=0'
      READ(5,*)IWRIT
      PRINT*,'COMPONENT: Z=1,RADIAL=,ANTIPODES=3'
      READ(5,*)ICOMP
      IF(ICOMP.EQ.1)EPS=1.
      IF(ICOMP.EQ.2)EPS=-1.
      IF(ICOMP.EQ.3)EPS=0.
        read(5,*)nchans
      read(5,33)(ichan(i),i=1,nchans)
33     format(a6)
      READ(5,*)NSTAT
c      print*,'nstat',nstat
      DO 35 K=1,NSTAT
      DO 34 J=1,5
34    IND(K,J)=0
35      READ(5,36,END=40)(INM(K,J),J=1,4),(IND(K,J),J=1,4)
36      FORMAT(2A2,2A1,4I2)
40      CONTINUE
c       print*,'dist2 entered'
      CALL DISTZ2(INM,IND,NSTAT)
       call ahopen(luin,inahfile,0,istat)
        rewind 11
      READ(11,5)ITITRE,EPLAT,EPLONG
      WRITE(6,5)ITITRE,EPLAT,EPLONG
        iter=0
20       continue
       call ahrewind(luin,istat)
c        print*,'after rewind'
c       print*,'istat',istat
       READ(11,2,END=900)IST,STLAT,STLON,DIST,AZ,iph
2      FORMAT(a6,1X,4F10.2,1X,I2,4(F7.2,2X))
       WRITE(6,1)ist
10     continue
       call rdahhead(luin,stcod,chan,stla,stlo,stel,jyr,jmo,
     *jda,jho,jmi,asec,kyr,lmo,lj,ih1,im,sec,dtsamp,ndata,ierr)
       la=kyr-1900
      if(ierr.ne.0)go to 20
c      write(6,104)stcod,chan,dtsamp,ndata,stla,stlo,jyr,jda,jho,jmi
104    format('station',a6,' channel:',a6,' dt:',f6.2,' np:',i6,
     *    /,'latitude:',f7.2,' longitude:',f8.2,4i4,f7.2)
c       print*,'ist stcod',ist,stcod
c        print*,'chan ichan',chan,ichan
c       print*,'jda jho jmi asec',jda,jho,jmi,asec
c        print*,'lj lh1 im sec',lj,lh1,im,sec
       if(ist(1:3).ne.stcod(1:3))go to 10
        do 144 ic=1,nchans
c	print*,'chan',chan,ichan(ic)
        if(chan(1:4).eq.ichan(ic)(1:4))go to 145
144      continue
        goto 10
145     continue
	print*,'channel',chan(1:4)
       call rdahdata(luin,x,ndata,ierr)
c        write(6,115)(x(i),i=1,200)
115      format(10e12.3)
        ih0=jho
       ih=ih1
c     if(jda.eq.(lj-1))ih=ih1+24
c       print*,'jda lj ih ih0',jda,lj,ih,ih0
      im0=jmi
      sl=asec
      dt=dtsamp
c	print*,'dt=',dt
	if(abs(dt-10.).lt.1.e-3)dt=10.
      call julien(jda,jmo,jyr,iy0,id0,1)
       l=ndata
       n=l
      U1(IPH)=U10(IPH)
      U2(IPH)=U20(IPH)
	print*,'u1 u2 iph',u1(iph),u2(iph),iph
      PRINT*,'jda,jMO,jyr,IY0,ID0',jda,jMO,jyr,IY0,ID0
C
C    CALCUL DE LA FENETRE EN TEMPS
C
37    NST=0
	print*,'l n',l,n
      DO 367 K=1,L
367    XW(K)=X(K)
	print*,'u1 u2 iph',u1(iph),u2(iph),iph
      T1=DIST/U1(IPH)
      T2=DIST/U2(IPH)
	if(iwrit.eq.0)go to 368
c	print*,'u1 u2 befor',iph, u1(iph),u2(iph)
c
c adjust group velocity window for purpose of tapering
c
	n1=t1/dt
	n2=t2/dt
	n=n2-n1+1
	nu=n1-n/2
	nv=n2+n/2
	U1(iph)=dist/(nu*dt)
	u2(iph)=dist/(nv*dt)
	t1=dist/u1(iph)
	t2=dist/u2(iph)
c	print*,'u1 u2 after',iph,u1(iph),u2(iph)
368	continue
      DELAY=SEC-SL+60.*FLOAT(IM-IM0)+3600.*FLOAT(IH-IH0)+
     *86400.*(lj-jda)
      DELAY=-DELAY
       WRITE(6,132)DELAY,SEC,sl,IS0,IM,IM0,IH,IH0
132       FORMAT(' DELAY SEC',3e12.3,5I5)
      N0=(T1-DELAY)/DT-NST
	print*,'n0,ndata',n0,ndata
	if(n0.ge.ndata)go to 10
       iter=iter+1
      N1=(T2-T1)/DT
      TS=N0*DT+DELAY+DT*NST
	print*,'t1 t2 n0 n1',t1,t2,n0,n1
c      print*,'dist ,t1 t2 u1(iph),u2(iph)',dist,t1,t2,u1(iph),u2(iph)
c      print*,'n0 n1 ts dt delay',n0,n1,ts,dt,delay
C     LX1=1024
      LX1=1024
C     NP=10
      NP=10
      MICO=DT-20.
      IF(MICO.EQ.0)THEN
                   LX1=256
                   NP=8
      ENDIF
19        CONTINUE
C
C   DETERMINATION DU NOMBRE DE POINTS POUR LA TF (PUISSANCE DE 2)
      NOUT=500
22      CONTINUE
      IF(N1.LE.LX1)GO TO 23
      LX1=2*LX1
      NP=NP+1
      NOUT=NOUT*2
      GO TO 22
23     CONTINUE
c      WRITE(6,333)LX1,NP,N0,N1,TS
333      FORMAT(' LX1,NP,N0,N1,TS',4I5,F9.2)
200     FORMAT(I5/(6E12.3))
C     WRITE(1,5)ITITRE
5         FORMAT(10A2,2F8.2)
4     FORMAT(6E12.5)
      IF(ICOMP.EQ.3)IPH=1+IPH
      WRITE(6,1)IST,STLAT,STLON,DIST,AZ,IPH,
     *DELAY,U1(IPH),U2(IPH),TS
C     DETERMINATION DES PERIODES A CONSERVER
C     SUR LA BASE DE TLF=LX1*DT
C
24      TLF=LX1*DT
      IF(TLF.GT.2000.)GO TO 25
      LX1=2*LX1
      NP=NP+1
      NOUT=NOUT*2
      GO TO 24
25     CONTINUE
      T(1)=99999.9
c	print*,'tlf',tlf
      DO 63 I=2,2048
63     T(I)=TLF/FLOAT(I-1)
      JJ=0
      DO 64 I=1,2048
      IF(T(I).GT.TT1)GO TO 64
      IF(T(I).LT.TT2)GO TO 65
      JJ=JJ+1
      TB(JJ)=T(I)
64     continue
65     CONTINUE
      NTOT=JJ
c      WRITE(6,105)NTOT
      IF(ITER.EQ.1)WRITE(9,105)NTOT
105      FORMAT(I5)
C     IF(ITER.EQ.1)WRITE(1,102)(TB(I),I=1,NTOT)
      IF(ITER.EQ.1)WRITE(9,102)(TB(I),I=1,NTOT)
      T(1)=0.0
      IUGE=1
1      FORMAT(a6,4F10.2,1X,I2,4(F7.2,2X))
503   DO 59 I=1,N1
      J=N0+I
      T(I+1)=T(I)+DT
CCC   WRITE(6,300)I,J,X(I),X(J)
300     FORMAT(' I J X',2I5,2E12.3)
59    YY(I)=XW(J)
C     GO TO 110
C     WRITE(6,220)(YY(I),I=1,20)
C     WRITE(6,220)(YY(I),I=N1-15,N1)
C      WRITE(6,100) N,N0,N1,LX1,NP,NOUT,DT
C     WRITE(6,200)N1,(X(I),I=1,N1)
      CALL DTRD(YY,N1,X,0,AMOY)
C     WRITE(1,220)(X(I),I=1,N1)
220    FORMAT(' x'/(6E12.3))
      PRINT*,'AMOY',AMOY
      PRINT*,'IUGE',IUGE
      IF(IUGE.EQ.0)GO TO 509
      IF(N0.EQ.0)GO TO 509
      DI0=ABS(X(1))
      DIF0=ABS(XW(N0)-AMOY)
C     PRINT*,'DI0',DI0,'DIF0',DIF0
C     PRINT*,'X(1)',X(1),'XW',XW(N0),'U1',U1(IPH)
      IF(DI0.GT.DIF0)THEN
                     N0=N0-1
                     N1=N1+1
                      U1(IPH)=DIST/(N0*DT+DELAY)
                         TS=N0*DT+DELAY+DT*NST
                     GO TO 503
                     ELSE
                DI0=ABS(X(1))
                DIF0=ABS(X(2))
C     PRINT*,'X1 X2',X(1),X(2),X(3),X(4),X(5)
              IF(DI0.GT.DIF0)THEN
              N0=N0+1
      N1=N1-1
      U1(IPH)=DIST/(N0*DT+DELAY)
                  TS=N0*DT+DELAY+DT*NST
              GO TO 503
               ELSE
                     IUGE=0
      ENDIF
                     ENDIF
C     PRINT*,'N0',N0
 509  DI1=ABS(X(N1))
         DIF=ABS(XW(N0+N1+1)-AMOY)
C     PRINT*,'DI1',DI1,'DIF,U1',U1(IPH),DIF
C     PRINT*,'X(N1) XW(N0N11)',X(N1),XW(N0+N1+1)
      IF(DI1.GT.DIF)THEN
                    N1=N1+1
                    U2(IPH)=DIST/((N0+N1)*DT+DELAY)
                    IF(U2(IPH).LT.2.5)GO TO 507
                    GO TO 503
                    ELSE
                    DI1=ABS(X(N1))
                    DIF=ABS(X(N1-1))
C     PRINT*,'XN XNM1',X(N1),X(N1-1),X(N1-2),X(N1-3)
                    IF(DI1.GT.DIF)THEN
                    N1=N1-1
                    U2(IPH)=DIST/((N0+N1)*DT+DELAY)
                    IF(U2(IPH).LT.2.5)GO TO 507
                    GO TO 509
                    ENDIF
      ENDIF
110   CONTINUE
507   DO 501 I=1,N1
      J=N0+I
501   X(I)=XW(J)
      IF(IWRIT.EQ.0)GO TO 511
	print*,n1,U1(iph),u2(iph)
      WRITE(9,1004)N1,U1(IPH),U2(IPH)
1004  FORMAT(I6,2F5.2)
      WRITE(9,102)(X(I),I=1,N1)
511   CONTINUE
C     CALL TAPER(X,N1,0.,0.5,X)
C     WRITE(6,102)(YY(I),I=1,N1)
      DO 500 IGNA=1,N1
500   YY(IGNA)=X(IGNA)
      DO 510 I=N1+1,LX1
      YY(I)=0.
510   CONTINUE
c	print*,'n1,lx1',n1,lx1
11      FORMAT(36A2)
103   FORMAT(6G12.5)
100   FORMAT(6I7,5G12.5)
C     WRITE(6,311)N1
311        FORMAT(' N1',I5)
      IF(IPLOT.NE.1)GO TO 341
C
C   CALCUL DE LA TRANSFORMEE DE FOURIER APRES ELIMINATION
C   DE LA DERIVE ET APODISATION
C
341      CONTINUE
c     BEGIN=0.25
c      begin=0.2
c      END=0.25
	begin=0.2
	end=0.2
      CALL DTRD(YY,N1,YY,0,AMIO)
      CALL TAPER(YY,N1,BEGIN,END,Y)
C      APODISATION FACON OUIZA
C     SOM=0.
C     DO 104 I=1,N1
C104    SOM=SOM+YY(I)
C     SOM=SOM/N1
C     DO 106 I=1,N1
C106     YY(I)=YY(I)-SOM
C     MI=(N1+1)/2
C     MT=2*MI
C     DO 114 I=1,N1
C     DI=I
C      DMI=MI
C114      Y(I)=YY(I)*(1.-(((DI-MI)**2)/(DMI**2))**2)
CC     PRINT*,'Y',(Y(I),I=1,N1)
      DO 504 KTI=N1+1,LX1
504   Y(KTI)=0.
      TS=N0*DT+DELAY+DT*NST
      ALX1=LX1
      TLF=ALX1*DT
      DO 51 I=1,LX1
      C(I)=CMPLX(Y(I), 0.0)
51    CONTINUE
      CALL COOL(NP,C,-1.0)
      NSP=NOUT
C     WRITE(6,444)NSP
444    FORMAT(' NSP',I5)
      AM=FLOAT(IPH)-1.
12       FORMAT(A6,i4,I3,F10.3,F5.1,I6,2F10.3,F5.0,2F8.2)
      DO 52 I=2, LX1
      Y(I)=CABS(C(I))*DT
      YY(I)=ATAN2(AIMAG(C(I)),REAL(C(I)))
52    CONTINUE
      T(1)=99999.9
      DO 53 I=2, LX1
      AI=I
      T(I)=TLF/(AI-1.0)
53    CONTINUE
C
C    SELECTION DES PERIODES CONSERVEES
C
      JJ=0
      DO 54 I=1,LX1
      IF(T(I).GT.TT1)GO TO 54
      IF(T(I).LT.TT2)GO TO 55
C     WRITE(6,446)T(I),TT1,TT2
446      FORMAT(' T TT1 TT2',3F10.3)
C     IF(TLF-10240.)57,56,57
57       DO 58 K=1,NTOT
      DIFF=ABS(T(I)-TB(K))
C     WRITE(1,447)K,T(I),TB(K),DIFF
      IF(DIFF.LT.1.E-3)GO TO 56
58       CONTINUE
447     FORMAT(' K T TB DIFF ',I5,2F10.3,E12.4)
      GO TO 54
56       CONTINUE
C     IF(I.LT.N1)GO TO 54
C     IF(I.GT.2*N1)GO TO 55
      JJ=JJ+1
      TT(JJ)=T(I)
      X(JJ)=Y(I)
      XX(JJ)=YY(I)
54      CONTINUE
55       CONTINUE
      NTOT=JJ
c      WRITE(6,445)NTOT
445     FORMAT(' NTOT',I5)
C     WRITE(1,105)NTOT
C      WRITE(1,102)(TT(I),I=1,NTOT)
C     WRITE(9,1004)N1,U1,U2
C1004  FORMAT(I6,2F5.2)
C     WRITE(9,102)(XW(I),I=1,N1)
c     WRITE(6,101)(TT(I),X(I),XX(I),I=1,NTOT)
C
C      CORRECTION INSTRUMENTALE ET DISPERSION GEOMETRIQUE
C      LE SPECTRE CORRIGE CORRESPOND A UN DEPLACEMENT EN CM
c
	jchn=0
	print*,'channel',chan(1:3)
c       if(chan(1:3).ne.'VLP'.and.chan(1:3).ne.'LP_')then
c	if(chan(1:3).ne.'VBB'.or.chan(1:3).ne.'IDA')go to 128
c	endif
        print*,'VLP,LP,VBB,IDA'
c       call getahpz(ds,a0,np,nz,pole,zero)
c	get response from ah header
	np=xnpah
	nz=xnzah
	ds=dsah
	a0=a0ah
	prod=ds*a0
	do i=1,np
	pole(i)=pz(1,i)
	enddo
	do i=1,nz
	zero(i)=pz(2,i)
	enddo
       print*,'ds a0 np nz',ds,a0,np,nz
       write(6,127)(pole(i),i=1,np)
        write(6,127)(zero(i),i=1,nz)
127     format(6e12.3) 
      DO 27 L=1,NTOT
      PER=TT(L)
      FR=1./PER
       om=fr*tpi
      FAINR=0
      GR=1.
       call digresp(om,nz,zero,np,pole,rea,aim)
c	convert from m to cm (assuming ahheader has displacement in meters)
       gr=rea*rea+aim*aim
	gr=sqrt(gr)*ds*a0*1.e-2
       fainr=atan2(aim,rea)
      FAINR=-FAINR
      FAIT0=-TPI*TS/PER
      FAIPS=-AM*PI/2.
      AC=ABS(SIN(DIST*PI/20000.0))
      AC=SQRT(AC)
      X(L)=X(L)*AC/GR
	print*,'per fait0 faips fainr eps',per,fait0,faips,fainr,eps
      CRFAI=XX(L)+FAIT0+FAIPS+FAINR-EPS*PI/4.
      CRFAI=CRFAI/TPI
      XX(L)=(CRFAI-AINT(CRFAI))*TPI
27        CONTINUE
C
C     ECRITURE SUR L'UNITE 9 : SPECTRES CORRIGES POUR L'INSTRUMENT ET
C     LA DISPERSION GEOMETRIQUE
C
C    ECRITURE DES DONNEES EN BINAIRE
      WRITE(9,12)IST,IY0,ID0,TS,DT,NTOT,DIST,AZ,AM
     *,STLAT,STLON
102      FORMAT(6(1X,E12.5))
      WRITE(9,102)(X(I),XX(I),I=1,NTOT)
101   FORMAT( 3(3G12.5,3X))
134       CONTINUE
C
C
C
c       print*,'going back to 20'
      GO TO 20
900       CONTINUE
       close(9)
      CLOSE(15)
      stop
      END
      FUNCTION ITA(LEOSTA)
      INTEGER*2 LEOSTA,IA(5)
      DATA IA/'Z','L','T','N','E'/
      ITA=1
      DO 10 I=1,5
      IF(LEOSTA.EQ.IA(I))GO TO 4
10    CONTINUE
      ITA=0
4     RETURN
      END
      SUBROUTINE JULIEN(IJ,IM,IA,IY,ID,ISI)
C
C     CONVERTI IJ,IM,IA -> IY,ID       POUR ISI=+1
C     CONVERTI IY,ID    -> IJ,IM,IA    POUR ISI=-1
C
C     EX:  IJ,IM,IA = 21 02 80    IY,ID = 1980 52
C
c      INTEGER*2 IJ,IM,IA,IY,ID,ISI
      DIMENSION IMO(12)
      DATA IMO/31,28,31,30,31,30,31,31,30,31,30,31/
      IF(ISI.EQ.-1)GO TO 10
      IF(IA.EQ.0)GO TO 20
      IMO(2)=28
      IF(MOD(IA,4).EQ.0)IMO(2)=29
      ID=0
      DO 1 I=1,12
      IF(I.EQ.IM)GO TO 2
1     ID=ID+IMO(I)
2     ID=ID+IJ
c      IY=1900+IA
       iy=ia
      RETURN
10    IF(IY.EQ.0)GO TO 20
      IMO(2)=28
      IF(MOD(IY,4).EQ.0)IMO(2)=29
      IDL=ID
      DO 3 I=1,12
      IF(IDL.LE.IMO(I))GO TO 4
3     IDL=IDL-IMO(I)
4     IJ=IDL
      IM=I
      IA=IY-1900
      RETURN
20    IJ=0
      IM=0
      IA=0
      IY=0
      ID=0
      RETURN
      END
      SUBROUTINE ITTA(LEOSTA,IST1,ITA,ITO)
C    IO=V MEANS THAT STATION RESPONSE IS IN VELOCITY, OTHERWISE ACCELERATION
      INTEGER*2 LEOSTA,IA(6),IST1,IO
      DATA IA/'Z','L','T','N','E','O'/
      DATA IO/'V'/
      ITA=1
      ITO=1
C     WRITE(1,555)LEOSTA
C    ITA=0    IDA
C    ITA=1    GEOSCOPE
C    ITA=2    SRO (RECHERCHE un Fichier SROPOLE)
555    FORMAT(' LEOSTA',1X,A1)
      DO 10 I=1,5
      IF(LEOSTA.EQ.IA(I))GO TO 4
10    CONTINUE
      IF(LEOSTA.NE.IA(6))GO TO 8
      ITA=2
      GO TO 4
8       CONTINUE
      ITA=0
4       CONTINUE
      IF(IST1.EQ.IO)GO TO 5
   20 CONTINUE
      ITO=0
    5 CONTINUE
C     WRITE(1,200)IST1,ITO,ITA
200     FORMAT(' IST1 ITO ITA',A1,2I5)
      CLOSE(12)
      RETURN
      END