c     FLAG=0 : FORWARD FFT CALCULATION
      implicit real(a-h,o-z)
      PARAMETER(NUM=2048,MEV=21,MFQ=17)
      REAL W(NUM),XR1(NUM),XR2(NUM),XI1(NUM),XI2(NUM),AN1(NUM),
     &     AN2(NUM),X1(2001),X2(2001),F(NUM),AN(num),ave(MFQ),ans1(num),
     &     ans(num),s(MFQ),azs(1),azc(1),rt(MFQ),tra(MFQ),rad(MFQ),
     &     r(MFQ),
     &     tra1(NUM),trar(NUM),trai(MFQ),t(MFQ),
     &     rad1(NUM),radr(NUM),radi(MFQ)
      complex XX1(NUM),XX2(NUM)
      character a*10,infileTOP*100,infileEND1*5,infileEND2*5,
     &          infileEND3*6,
     &          INFILE1(MEV)*100,INFILE2(MEV)*100,INFILE3(MEV)*100,
     &          OTFILE1(MEV)*100,OTFILE2(MEV)*100,
     &          OTFILE3*100,OTFILE4*100,infileDATE*10,
     &          otfileTOP*100,otfileEND1*5,otfileEND2*5,
     &          OTFILE5*100,OTFILE6*100
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      NEV=MEV

      PI=4.0*ATAN(1.0)

      N=2048

      DT=0.005

      DF=1./(N*DT)

      a='0123456789'

      open (92,file='../result/e16.inp')

      do 913 ID2=1,16
c      write(*,*)'980114   [ 1 ]'
c      write(*,*)'980116   [ 2 ]'
c      write(*,*)'980516   [ 3 ]'
c      write(*,*)'980829   [ 4 ]'
c      write(*,*)'981108   [ 5 ]'
c      write(*,*)'981128   [ 6 ]'
c      write(*,*)'981203   [ 7 ]'
c      write(*,*)'99042518 [ 8 ]'
c      write(*,*)'990715   [ 12 ]'
c      write(*,*)'990809   [ 9 ]'
c      write(*,*)'990811   [ 10 ]'
c      write(*,*)'990913   [ 11 ]'
c      write(*,*)'000410   [ 13 ]'
c      write(*,*)'000603   [ 14 ]'
c      write(*,*)'000827   [ 15 ]'
c      write(*,*)'000929   [ 16 ]'

      if (ID2.eq.1) then
         infileDATE='980114'
         OTFILE5='../result/980114/data1.out'
         OTFILE6='../result/980114/ratio1.out'
      elseif (ID2.eq.2) then
         infileDATE='980116'
         OTFILE5='../result/980116/data2.out'
         OTFILE6='../result/980116/ratio2.out'
      elseif (ID2.eq.3) then
         infileDATE='980516'
         OTFILE5='../result/980516/data3.out'
         OTFILE6='../result/980516/ratio3.out'
      elseif (ID2.eq.4) then
         infileDATE='980829'
         OTFILE5='../result/980829/data4.out'
         OTFILE6='../result/980829/ratio4.out'
      elseif (ID2.eq.5) then
         infileDATE='981108'
         OTFILE5='../result/981108/data5.out'
         OTFILE6='../result/981108/ratio5.out'
      elseif (ID2.eq.6) then
         infileDATE='981128'
         OTFILE5='../result/981128/data6.out'
         OTFILE6='../result/981128/ratio6.out'
      elseif (ID2.eq.7) then
         infileDATE='981203'
         OTFILE5='../result/981203/data7.out'
         OTFILE6='../result/981203/ratio7.out'
      elseif (ID2.eq.8) then
         infileDATE='99042518'
         OTFILE5='../result/99042518/data8.out'
         OTFILE6='../result/99042518/ratio8.out'
      elseif (ID2.eq.12) then
         infileDATE='990715'
         OTFILE5='../result/990715/data17.out'
         OTFILE6='../result/990715/ratio17.out'
      elseif (ID2.eq.9) then
         infileDATE='990809'
         OTFILE5='../result/990809/data9.out'
         OTFILE6='../result/990809/ratio9.out'
      elseif (ID2.eq.10) then
         infileDATE='990811'
         OTFILE5='../result/990811/data10.out'
         OTFILE6='../result/990811/ratio10.out'
      elseif (ID2.eq.11) then
         infileDATE='990913'
         OTFILE5='../result/990913/data11.out'
         OTFILE6='../result/990913/ratio11.out'
      elseif (ID2.eq.13) then
         infileDATE='000410'
         OTFILE5='../result/000410/data18.out'
         OTFILE6='../result/000410/ratio18.out'
      elseif (ID2.eq.14) then
         infileDATE='000603'
         OTFILE5='../result/000603/data19.out'
         OTFILE6='../result/000603/ratio19.out'
      elseif (ID2.eq.15) then
         infileDATE='000827'
         OTFILE5='../result/000827/data20.out'
         OTFILE6='../result/000827/ratio20.out'
      elseif (ID2.eq.16) then
         infileDATE='000929'
         OTFILE5='../result/000929/data21.out'
         OTFILE6='../result/000929/ratio21.out'
      endif

      open (93,file=OTFILE5)
      open (94,file=OTFILE6)


      call CountNumberOfChracters2(infileDATE,1,ie00)

      do 912 ID=1,18
c      write(*,*)'aoba     [ 13 ]'
c      write(*,*)'asahi    [ 8 ]'
c      write(*,*)'hodogaya [ 7 ]'
c      write(*,*)'isogo    [ 9 ]'
c      write(*,*)'izumi    [ 17 ]'
c      write(*,*)'kanagawa [ 2 ]'
c      write(*,*)'kanazawa [ 10 ]'
c      write(*,*)'kouhoku  [ 11]'
c      write(*,*)'kounan   [ 6 ]'
c      write(*,*)'midori   [ 12 ]'
c      write(*,*)'minami   [ 5 ]'
c      write(*,*)'naka     [ 4 ]'
c      write(*,*)'nishi    [ 3 ]'
c      write(*,*)'sakae    [ 16 ]'
c      write(*,*)'seya     [ 18 ]'
c      write(*,*)'totuka   [ 15 ]'
c      write(*,*)'tuduki   [ 14 ]'
c      write(*,*)'turumi   [ 1 ]'
c      read(*,*)ID

      if (ID.eq.13) then
         infileTOP='../'//infileDATE(1:ie00)//'/aoba/ab'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/aoba/ab'
      elseif (ID.eq.8) then
         infileTOP='../'//infileDATE(1:ie00)//'/asahi/as'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/asahi/as'
      elseif (ID.eq.7) then
         infileTOP='../'//infileDATE(1:ie00)//'/hodogaya/hd'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/hodogaya/hd'
      elseif (ID.eq.9) then
         infileTOP='../'//infileDATE(1:ie00)//'/isogo/is'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/isogo/is'
      elseif (ID.eq.17) then
         infileTOP='../'//infileDATE(1:ie00)//'/izumi/iz'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/izumi/iz'
      elseif (ID.eq.2) then
         infileTOP='../'//infileDATE(1:ie00)//'/kanagawa/kg'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/kanagawa/kg'
      elseif (ID.eq.10) then
         infileTOP='../'//infileDATE(1:ie00)//'/kanazawa/kz'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/kanazawa/kz'
      elseif (ID.eq.11) then
         infileTOP='../'//infileDATE(1:ie00)//'/kouhoku/kh'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/kouhoku/kh'
      elseif (ID.eq.6) then
         infileTOP='../'//infileDATE(1:ie00)//'/kounan/kn'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/kounan/kn'
      elseif (ID.eq.12) then
         infileTOP='../'//infileDATE(1:ie00)//'/midori/md'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/midori/md'
      elseif (ID.eq.5) then
         infileTOP='../'//infileDATE(1:ie00)//'/minami/mi'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/minami/mi'
      elseif (ID.eq.4) then
         infileTOP='../'//infileDATE(1:ie00)//'/naka/na'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/naka/na'
      elseif (ID.eq.3) then
         infileTOP='../'//infileDATE(1:ie00)//'/nishi/ns'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/nishi/ns'
      elseif (ID.eq.16) then
         infileTOP='../'//infileDATE(1:ie00)//'/sakae/sk'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/sakae/sk'
      elseif (ID.eq.18) then
         infileTOP='../'//infileDATE(1:ie00)//'/seya/sy'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/seya/sy'
      elseif (ID.eq.15) then
         infileTOP='../'//infileDATE(1:ie00)//'/totuka/tt'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/totuka/tt'
      elseif (ID.eq.14) then
         infileTOP='../'//infileDATE(1:ie00)//'/tuduki/td'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/tuduki/td'
      elseif (ID.eq.1) then
         infileTOP='../'//infileDATE(1:ie00)//'/turumi/tr'
         otfileTOP='../result/'//infileDATE(1:ie00)//'/turumi/tr'
      endif


      call CountNumberOfChracters2(infileTOP,1,ie01)
      call CountNumberOfChracters2(otfileTOP,1,ie02)


      infileEND1='e.inp'
      infileEND2='n.inp'
      infileEND3='az.inp'
      otfileEND1='.out1'
      otfileEND2='.rat1'

      OTFILE3=otfileTOP(1:ie02)//'.out1'
      OTFILE4=otfileTOP(1:ie02)//'.rat1'

      open(67,file=OTFILE3)
      open(68,file=OTFILE4)

      do 586 iEV=1,NEV


cccc ループから外すもの（対象外） cccc
         if (iEV.ge.22) then
            goto 586
         endif
cccccccccccccccccccccccccccccccccccccc

         if (iEV.lt.10) then
      INFILE1(iEV)=infileTOP(1:ie01)//'0'//a(iEV+1:iEV+1)//infileEND1
      INFILE2(iEV)=infileTOP(1:ie01)//'0'//a(iEV+1:iEV+1)//infileEND2
      INFILE3(iEV)=infileTOP(1:ie01)//'0'//a(iEV+1:iEV+1)//infileEND3
      OTFILE1(iEV)=otfileTOP(1:ie02)//'0'//a(iEV+1:iEV+1)//otfileEND1
      OTFILE2(iEV)=otfileTOP(1:ie02)//'0'//a(iEV+1:iEV+1)//otfileEND2
         else if (iEV.ge.10.and.iEV.lt.20) then
      INFILE1(iEV)=infileTOP(1:ie01)//'1'//a(iEV-9:iEV-9)//infileEND1
      INFILE2(iEV)=infileTOP(1:ie01)//'1'//a(iEV-9:iEV-9)//infileEND2
      INFILE3(iEV)=infileTOP(1:ie01)//'1'//a(iEV-9:iEV-9)//infileEND3
      OTFILE1(iEV)=otfileTOP(1:ie02)//'1'//a(iEV-9:iEV-9)//otfileEND1
      OTFILE2(iEV)=otfileTOP(1:ie02)//'1'//a(iEV-9:iEV-9)//otfileEND2
         else if (iEV.ge.20.and.iEV.lt.30) then
      INFILE1(iEV)=infileTOP(1:ie01)//'2'//a(iEV-19:iEV-19)//infileEND1
      INFILE2(iEV)=infileTOP(1:ie01)//'2'//a(iEV-19:iEV-19)//infileEND2
      INFILE3(iEV)=infileTOP(1:ie01)//'2'//a(iEV-19:iEV-19)//infileEND3
      OTFILE1(iEV)=otfileTOP(1:ie02)//'2'//a(iEV-19:iEV-19)//otfileEND1
      OTFILE2(iEV)=otfileTOP(1:ie02)//'2'//a(iEV-19:iEV-19)//otfileEND2
         endif
         call CountNumberOfChracters(INFILE1,iEV,1,iend1,MEV)
         call CountNumberOfChracters(INFILE2,iEV,1,iend2,MEV)
         call CountNumberOfChracters(INFILE3,iEV,1,iend3,MEV)
         call CountNumberOfChracters(OTFILE1,iEV,1,iend4,MEV)
         call CountNumberOfChracters(OTFILE2,iEV,1,iend5,MEV)

         open(11,file=INFILE1(iEV)(1:iend1),status='old',err=586)
         open(12,file=INFILE2(iEV)(1:iend2),status='old',err=586)
         open(21,file=INFILE3(iEV)(1:iend3),status='old',err=586)
         open(30,file=OTFILE1(iEV)(1:iend4))
         open(31,file=OTFILE2(iEV)(1:iend5))

c      write(*,*)OTFILE1(iEV)(3:40)
c      write(*,*)OTFILE2(iEV)(3:40)


c      OPEN(11, file='INDATA1',status='old')
c      OPEN(12, file='INDATA2',status='old')
c      open(21, file='indata3',status='old')

      DO I=1,30
         READ(11,*)
      ENDDO
      READ(11,*)(X1(L),L=1,2001)
      close(11)

      DO J=1,30
         READ(12,*)
      ENDDO
      READ(12,*)(X2(J),J=1,2001)
      close(12)

      read(21,*)azs(1)
      read(21,*)azc(1)
      close(21)

      DO 1 K=1,2001
         XR1(K)=X1(K)
         XR2(K)=X2(K)
1     CONTINUE
      DO 2 L=2002,2048
         XR1(L)=0.0
         XR2(L)=0.0
2     CONTINUE

      DO 10 I=1,NUM
         W(I)=2.*PI*FLOAT(I)/(FLOAT(N)*DT)
10    continue

c     FORWARD FFT

      CALL FFT1(XR1,XI1,N,ITER,0)
      CALL FFT1(XR2,XI2,N,ITER,0)

c      do 571 i=1,NUM
c         XX1(i)=cmplx(XR1(i),XI1(i))
c         XX2(i)=cmplx(XR2(i),XI2(i))
c571   continue
c      call nlogn(11,XX1,-1.)
c      call nlogn(11,XX2,-1.)
c      do 572 i=1,NUM
c         XX1(i)=XX1(i)*DT
c         XX2(i)=XX2(i)*DT
c         XR1(i)=real(XX1(i))
c         XI1(i)=imag(XX1(i))
c         XR2(i)=real(XX2(i))
c         XI2(i)=imag(XX2(i))
c572   continue



      DO 21 I=1,NUM
         AN1(I)=SQRT(XR1(I)**2.+XI1(I)**2.)
21    CONTINUE
      DO 22 J=1,NUM
         AN2(J)=SQRT(XR2(J)**2.+XI2(J)**2.)
22    CONTINUE
      DO 23 K=1,NUM
         ans(K)=SQRT(AN1(K)**2.+AN2(K)**2.) 
         F(K)=K*DF
23    CONTINUE
c      CALL SWIN(2048,ans,G1,2048,100,DF,0.390625)

      do k=1,num
         an(k)=ans(k)
      enddo


ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       do 100 k=1,17
c            s(k)=0
c 100        continue
c       do 101 i=1,7
c             s(1)=s(1)+an(i+1)
c             s(2)=s(2)+an(i+4) 
c             s(3)=s(3)+an(i+7)
c             s(4)=s(4)+an(i+10)
c             s(5)=s(5)+an(i+13)
c             s(6)=s(6)+an(i+16)
c 101         continue
c
c       do 102 j=1,11
c              s(7)=s(7)+an(j+19)
c              s(8)=s(8)+an(j+24)
c              s(9)=s(9)+an(j+29)
c              s(10)=s(10)+an(i+34)
c              s(11)=s(11)+an(i+39)
c              s(12)=s(12)+an(i+44)
c 102      continue
c
c       do 106 j=1,21
c              s(13)=s(13)+an(j+49)
c              s(14)=s(14)+an(j+59)
c              s(15)=s(15)+an(j+69)
c              s(16)=s(16)+an(i+79)
c              s(17)=s(17)+an(i+89)
c 106          continue
c
c         do l=1,6
c             ave(l)=(s(l))/7
c           enddo  
c         do l=7,12
c             ave(l)=(s(l))/11
c             enddo
c         do l=13,17
c              ave(l)=(s(l))/21
c              enddo
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
      do 100 k=1,17
         s(k)=0
100   continue

      do 101 j=1,3
         s(1)=s(1)+an(j+3)
         s(2)=s(2)+an(j+6)
         s(3)=s(3)+an(j+8)
101   continue

      do 102 j=1,5
         s(4)=s(4)+an(j+10)
         s(5)=s(5)+an(j+12)
         s(6)=s(6)+an(j+17)
102   continue

      do 103 j=1,7
         s(7)=s(7)+an(j+21) 
         s(8)=s(8)+an(j+26)
103   continue

      do 104 i=1,9
         s(9)=s(9)+an(j+30)
         s(10)=s(10)+an(j+35)
104   continue

      do 183 i=1,11
         s(11)=s(11)+an(j+39)
         s(12)=s(12)+an(j+44)
183   continue

      do 184 i=1,13
         s(13)=s(13)+an(i+53)
184   continue

      do 185 i=1,15
         s(14)=s(14)+an(i+62)
185   continue

      do 186 i=1,17
         s(15)=s(15)+an(i+71)
186   continue

      do 187 i=1,19
         s(16)=s(16)+an(i+80)
187   continue

      do 188 i=1,21
         s(17)=s(17)+an(i+89)
188   continue

      ave(1)=s(1)/3
      ave(2)=s(2)/3
      ave(3)=s(3)/3
      ave(4)=s(4)/5
      ave(5)=s(5)/5
      ave(6)=s(6)/5
      ave(7)=s(7)/7
      ave(8)=s(8)/7
      ave(9)=s(9)/9
      ave(10)=s(10)/9
      ave(11)=s(11)/11
      ave(12)=s(12)/11
      ave(13)=s(13)/13
      ave(14)=s(14)/15
      ave(15)=s(15)/17
      ave(16)=s(16)/19
      ave(17)=s(17)/21
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      do 101 j=1,5
c         s(1)=s(1)+an(j+2)
c         s(2)=s(2)+an(j+4)
c101   continue
c
c      do 102 j=1,7
c         s(3)=s(3)+an(j+6)
c         s(4)=s(4)+an(j+9)
c         s(5)=s(5)+an(j+12)
c         s(6)=s(6)+an(j+16)
c102   continue
c
c      do 103 j=1,11
c         s(7)=s(7)+an(j+19) 
c         s(8)=s(8)+an(j+25)
c         s(9)=s(9)+an(j+30)
c         s(10)=s(10)+an(j+35)
c         s(11)=s(11)+an(j+40)
c         s(12)=s(12)+an(j+45)
c103   continue
c
c      do 104 i=1,21
c         s(13)=s(13)+an(i+50)
c         s(14)=s(14)+an(i+60)
c         s(15)=s(15)+an(i+70)
c         s(16)=s(16)+an(i+80)
c         s(17)=s(17)+an(i+90)
c104   continue
c
c      ave(1)=s(1)/5
c      ave(2)=s(2)/5
c      ave(3)=s(3)/7
c      ave(4)=s(4)/7
c      ave(5)=s(5)/7
c      ave(6)=s(6)/7
c      ave(7)=s(7)/11
c      ave(8)=s(8)/11
c      ave(9)=s(9)/11
c      ave(10)=s(10)/11
c      ave(11)=s(11)/11
c      ave(12)=s(12)/11
c      ave(13)=s(13)/21
c      ave(14)=s(14)/21
c      ave(15)=s(15)/21
c      ave(16)=s(16)/21
c      ave(17)=s(17)/21
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      do 105 j=1,11
c         s(11)=s(11)+an(j+39)
c         s(12)=s(12)+an(j+44)
c105   continue
c
c      do 106 j=1,13 
c         s(13)=s(13)+an(j+53)
c106   continue
c
c      do 107 i=1,15
c         s(14)=s(14)+an(j+62)
c107   continue
c
c      do 108 j=1,17
c         s(15)=s(15)+an(j+71)
c108   continue
c
c      do  109 j=1,19
c         s(16)=s(16)+an(80)
c109   continue
c
c      do 110 i=1,21
c         s(17)=s(17)+an(j+89)
c110   continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccc



      do i=1,NUM
         trar(i)=-AN1(i)*azc(1)+AN2(i)*azs(1)
         radr(i)=AN1(i)*azs(1)+AN2(i)*azc(1)
      enddo

c      CALL FFT1(trar,trai,N,ITER,0)
c      CALL FFT1(radr,radi,N,ITER,0)

      DO 31 I=1,NUM
         tra1(I)=SQRT(trar(I)**2.+trai(I)**2.)
31    CONTINUE
      DO 32 J=1,NUM
         rad1(J)=SQRT(radr(J)**2.+radi(J)**2.)
32    CONTINUE

      do 1100 k=1,17
         t(k)=0
         r(k)=0
1100  continue

c   radial component
      do 111 i=1,7
         r(1)=r(1)+rad1(i+1)
         r(2)=r(2)+rad1(i+4) 
         r(3)=r(3)+rad1(i+7)
         r(4)=r(4)+rad1(i+10)
         r(5)=r(5)+rad1(i+13)
         r(6)=r(6)+rad1(i+16)
111   continue

      do 112 j=1,11
         r(7)=r(7)+rad1(j+19)
         r(8)=r(8)+rad1(j+24)
         r(9)=r(9)+rad1(j+29)
         r(10)=r(10)+rad1(i+34)
         r(11)=r(11)+rad1(i+39)
         r(12)=r(12)+rad1(i+44)
112   continue

      do 116 j=1,21
         r(13)=r(13)+rad1(j+49)
         r(14)=r(14)+rad1(j+59)
         r(15)=r(15)+rad1(j+69)
         r(16)=r(16)+rad1(i+79)
         r(17)=r(17)+rad1(i+89)
116   continue

      do l=1,6
         rad(l)=(r(l))/7
      enddo

      do l=7,12
         rad(l)=(r(l))/11
      enddo

      do l=13,17
         rad(l)=(r(l))/21
      enddo

c  transverse component
      do 121 i=1,7
         t(1)=t(1)+tra1(i+1)
         t(2)=t(2)+tra1(i+4) 
         t(3)=t(3)+tra1(i+7)
         t(4)=t(4)+tra1(i+10)
         t(5)=t(5)+tra1(i+13)
         t(6)=t(6)+tra1(i+16)
121   continue

      do 122 j=1,11
         t(7)=t(7)+tra1(j+19)
         t(8)=t(8)+tra1(j+24)
         t(9)=t(9)+tra1(j+29)
         t(10)=t(10)+tra1(i+34)
         t(11)=t(11)+tra1(i+39)
         t(12)=t(12)+tra1(i+44)
122   continue

      do 126 j=1,21
         t(13)=t(13)+tra1(j+49)
         t(14)=t(14)+tra1(j+59)
         t(15)=t(15)+tra1(j+69)
         t(16)=t(16)+tra1(i+79)
         t(17)=t(17)+tra1(i+89)
126   continue

      do l=1,6
         tra(l)=(t(l))/7
      enddo

      do l=7,12
         tra(l)=(t(l))/11
      enddo

      do l=13,17
         tra(l)=(t(l))/21
      enddo

      do 200 k=1,17
         rt(k)=rad(k)/tra(k)
200   continue

c         write(6,*)tra1(1),rad1(1),rad(1),rt(1)
c         write(6,*)(rt(i),i=1,17)
c         write(6,*)(ave(i),i=1,17)

c      OPEN(30,file='OUTDATA')
c      open(31,file='outdata2')

      WRITE(30,1001)(ave(I),I=1,17)
      WRITE(67,1001)(ave(I),I=1,17)
      WRITE(93,1001)(ave(I),I=1,17)
      WRITE(92,1001)(ave(I),I=1,17)
      close(30)
      write(31,1001)(rt(l),l=1,17)
      write(68,1001)(rt(l),l=1,17)
      write(94,1001)(rt(l),l=1,17)
      close(31)
c100   format(17f8.5)
c      write(31,*)ave(1)
c      STAT = CODE(1:4) // '.OUT'
1000  FORMAT(105F8.4)
1001  FORMAT(17F10.5)


586   continue

      close(67)
      close(68)


912   continue

      close(93)
      close(94)

      write(*,*)infileDATE


913   continue

      close(92)


      STOP
      END

c******************************************************************************
c     1-DIMENSIONAL FFT BY COOLEY-TUKEY ALGORITHM
c******************************************************************************
      SUBROUTINE FFT1(XR,XI,N,ITER,FLAG)

      INTEGER FLAG,XP,XP2
      REAL XR(N),XI(N)

      IF (N.LT.2) THEN
        RETURN
      ENDIF
      IF (ITER.LE.0) THEN
        ITER=0
        I=N
   10   I=I/2
        IF (I.NE.0) THEN
          ITER=ITER+1
          GOTO 10
        ENDIF
      ENDIF
      J=1
      DO 20 I=1,ITER
        J=2*J
   20 CONTINUE
      IF (N.NE.J) THEN
        RETURN
      ENDIF
      IF (FLAG.EQ.0) THEN
        SGN=1.0
      ELSE
        SGN=-1.0
      ENDIF
      XP2=N
      DO 50 IT=1,ITER-1
        XP=XP2
        XP2=XP/2
        W=4.0*ATAN(1.0)/XP2
        DO 40 K=1,XP2
          ARG=(K-1)*W
          WR=COS(ARG)
          WI=SGN*SIN(ARG)
          I=K-XP
          DO 30 J=XP,N,XP
           J1=J+I
           J2=J1+XP2
           DR1=XR(J1)
           DR2=XR(J2)
           DI1=XI(J1)
           DI2=XI(J2)
           TR=DR1-DR2
           TI=DI1-DI2
           XR(J1)=DR1+DR2
           XI(J1)=DI1+DI2
           XR(J2)=TR*WR-TI*WI
           XI(J2)=TI*WR+TR*WI
   30     CONTINUE
   40   CONTINUE
   50 CONTINUE
      J1=N/2
      J2=N-1
      J=1
      DO 70 I=1,J2
        IF (I.LT.J) THEN
          TR=XR(J)
          TI=XI(J)
          XR(J)=XR(I)
          XI(J)=XI(I)
          XR(I)=TR
          XI(I)=TI
        ENDIF
        K=J1
   60    IF (K.LT.J) THEN
          J=J-K
          K=K/2
          GOTO 60
         ENDIF
         J=J+K
   70 CONTINUE
      IF (FLAG.EQ.0) THEN
        W=N
        DO 80 I=1,N
          XR(I)=XR(I)*10/W
          XI(I)=XI(I)*10/W
   80   CONTINUE
      ENDIF
      RETURN
      END
C************************************************
C  SUBROUTINE SMOOTHING
C************************************************
      SUBROUTINE SWIN(NFOLD,F,G,ND,IND,DF,BAND)
	DIMENSION F(ND),G(ND),W(101),G1(4497),G2(4497)
	T=1./DF
	IF(BAND.EQ.0.) GO TO 120
	UDF=1.854305/BAND*DF
	IF(UDF.GT.0.5) GO TO 230
	LMAX=IFIX(2./UDF)+1
C	IF(LMAX.GT.0.5) GO TO 240

	W(1)=0.75*UDF
	DO 110 L=2,LMAX
	DIF=1.570796*REAL(L-1)*UDF
	W(L)=W(1)*(SIN(DIF)/DIF)**4
110	CONTINUE
	
120	IF(IND.NE.100) GO TO 140
	G(1)=F(1)**2/T
	DO 130 K=2,NFOLD-1
	G(K)=2.*F(K)**2/T
130   CONTINUE
	G(NFOLD)=F(NFOLD)**2/T

140	IF(BAND.EQ.0.) GO TO 210
	LL=LMAX*2-1
	LN=LL-1+NFOLD
	LT=(LL-1)*2+NFOLD
	LE=LT-LMAX+1
	DO 150 K=1,LT
	G1(K)=0.
150	CONTINUE
	DO 160 K=1,NFOLD
	G1(LL-1+K)=G(K)
160   CONTINUE
      DO 180 K=LMAX,LE
	S=W(1)*G1(K)
	DO 170 L=2,LMAX
	S=S+W(L)*(G1(K-L+1)+G1(K+L-1))
170   CONTINUE
	G2(K)=S
180	CONTINUE
      DO 190 L=2,LMAX
	G2(LL+L-1)=G2(LL+L-1)+G2(LL-L+1)
	G2(LN-L+1)=G2(LN-L+1)+G2(LN+L+1)
190	CONTINUE
	DO 200 K=1,NFOLD
	G(K)=G2(LL-1+K)
200   CONTINUE
C
210	F(1)=SQRT(G(1)*T)
	DO 220 K=2,NFOLD-1
	F(K)=SQRT(G(K)*T/2.)
220	CONTINUE
	F(NFOLD)=SQRT(G(NFOLD)*T)
	RETURN
	
230   WRITE(6,501)
      STOP
240   WRITE(6,502)

501   FORMAT('BANDWIDTH IS TOO NARROW')
502   FORMAT('BANDWIDTH IS TOO WIDE')
      END

c******************************************************************************
c     nlogn
c******************************************************************************

      subroutine nlogn(n,x,sign)                                        
      complex   x,wk,hold,q                                            
      dimension m(15),x(2)                                              
      lx=2**n                                                           
      do 1 i=1,n                                                        
1     m(i)=2**(n-i)                                                     
      do 4 l=1,n                                                        
      nblock=2**(l-1)                                                   
      lblock=lx/nblock                                                  
      lbhalf=lblock/2                                                   
      k=0                                                               
      do 4 iblock=1,nblock                                              
      fk=k                                                              
      flx=lx                                                            
      v=sign*6.2831853e0*fk/flx                                         
      wk=cmplx(cos(v),sin(v))                                           
      istart=lblock*(iblock-1)                                          
      do 2 i=1,lbhalf                                                   
      j=istart+i                                                        
      jh=j+lbhalf                                                       
      q=x(jh)*wk                                                        
      x(jh)=x(j)-q                                                      
      x(j)=x(j)+q                                                       
2     continue                                                          
      do 3 i=2,n                                                        
      ii=i                                                              
      if(k.lt.m(i)) go to 4                                             
3     k=k-m(i)                                                          
4     k=k+m(ii)                                                         
      k=0                                                               
      do 7 j=1,lx                                                       
      if(k.lt.j) go to 5                                                
      hold=x(j)                                                         
      x(j)=x(k+1)                                                       
      x(k+1)=hold                                                       
5     do 6 i=1,n                                                        
      ii=i                                                              
      if(k.lt.m(i)) go to 7                                             
6     k=k-m(i)                                                          
7     k=k+m(ii)                                                         
      if(sign.lt.0.0) return                                            
      do 8 i=1,lx                                                       
8     x(i)=x(i)/flx                                                     
      return
      end




C***********************************************************************
      subroutine CountNumberOfChracters(a,n,istart,iend,MEV)
C***********************************************************************

      implicit real*8(a-h,o-z)

      character a(MEV)*100

      do 100 i=istart,100
         if (a(n)(i:i).eq.',') then
            go to 1
         else if (a(n)(i:i).eq.'') then
            go to 1
         end if
100   continue
1     iend=i-1

      return
      end

C***********************************************************************
      subroutine CountNumberOfChracters2(a,istart,iend)
C***********************************************************************

      implicit real*8(a-h,o-z)

      character a*100

      do 100 i=istart,100
         if (a(i:i).eq.',') then
            go to 1
         else if (a(i:i).eq.'') then
            go to 1
         end if
100   continue
1     iend=i-1

      return
      end




