c--gnoise.f
c--orders a number of files according to noise content
c---------------------------------------------------------------------
      program GNOISE
      parameter(ndim=180000)
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      dimension x(ndim),y(ndim),f(400),rm(400),
     &          a(10,10),b(10)
      character*60 title
      character*40 fname1,fname2,fname(400)
      character*20 fmt,dunit,tunit
      character*5 fday
      character*3 fext
      character*1 ans
      common // x,y
c
c--initialise
      ideg = 9
      maxp = ideg+1
c
c--files opened :
c  (input)
c       1 = data (residual file)
c  (output)
c       2 = data with polynomial subtracted
c       3 = ordered list of files
c
c--choose which files to be processed
      write(*,'(/a)') ' Program GNOISE'
      write(*,'(a)')  ' -------------------------------------------'
      write(*,'(a$)') ' Input filename extension (e.g. plm) --> '
      read(*,*) fext
      write(*,'(a$)') ' Enter first day to be processed     --> '
      read(*,*) iday1
      write(*,'(a$)') ' ..... last day ................     --> '
      read(*,*) iday2
c--if polynomial not already removed, ask about output data
      if(fext.ne.'dmp') then
        write(*,'(a$)') ' Write data minus polynomial [y/n]   --> '
        read(*,'(a1)') ans
        zwpoly = .false.
        if(ans.eq.'y') zwpoly=.true.
      endif
c
c--try to find specified files for the days specified
      kk = 0
      kday = iday1
c
c--convert kday to a character string 'fday '
   5  call intch(kday,fday)
c--then transfer fday to the input filename
      call fnday(fname1,'g',fday,fext)
c--check to see if input file exists
      inquire(file = fname1,exist =zexist)
      if(zexist) goto 8
      write(*,'(a,i5)') ' No input file found for day ',kday
      goto 10
c
c--copy the filename to the array
   8  kk = kk+1
      fname(kk) = fname1
c--make sure first good file begins the sequence
      if(kk.eq.1) iday1 = kday
c--try another day
  10  kday = kday+1
      if(kday.le.iday2) goto 5
c--there are kmax days to be processed
      kmax = kk
      k = 1
c
c--select a day from the list
  15  fname1 = fname(k)
c--give filename extension for output file and extract fday
      fname2 = fname1
      call fnext2(fname2,fday,'dmp')
      call chint(fday,kday)
c
c--read the data
      write(*,*) '  reading: ',fname(k)
      open(1,file=fname1)
      read(1,'(a)') title
      read(1,'(i10,g12.5)') np,delt
      read(1,'(3a20)') fmt,dunit,tunit
      read(1,fmt) (y(i),i=1,np)
      close(1)
c
c--if this is alread a 'dmp' file, then go directly to the statistics
      if(fext.eq.'dmp') goto 65
c
c--find and subtract 9th degree polynomial using pfit subroutine
c      write(*,*) '  fitting polynomial'
      do 20 i=1,np
        x(i) = float(i)
  20  continue
      call pfit(x,y,1,np,maxp,a,b)
c      write(*,'(5f12.5)') (b(i),i=1,maxp)
      do 60 ii=1,np
        sum = b(maxp)
        do 40 jj=maxp-1,1,-1
          sum = x(ii)*sum+b(jj)
  40    continue
        y(ii) = y(ii)-sum
  60  continue
c
c--write the data minus polynomial
      if(zwpoly) then
        write(*,*) '  writing: ',fname2
        open(2,file=fname2)
        write(2,'(a)') ' data minus polynomial'
        write(2,'(i10,g12.5)') np,delt
        write(2,'(3a20)') fmt,dunit,tunit
        write(2,fmt) (y(i),i=1,np)
        close(2)
      endif
c
c--remove mean
  65  call meanv(y,np,ymean,1)
c
c--find the rms deviation
      rms = 0.d0
      do 70 i=1,np
        rms = rms+(y(i)-ymean)**2
  70  continue
      rms = dsqrt(rms/float(np))
      f(k) =  float(kday)
      rm(k) = rms
      write(*,'(a8,i8,a5,f10.5)') '   day: ',kday,' rms ',rms
c
c--get another file
      k = k+1
      if(k.le.kmax) goto 15
c
c--sort the days according to rms
      call sort(rm,f,kmax)
c--write final output
      open(3,file='gnoise.out')
      write(3,'(a)') '     day      rms noise '
      write(3,'(a)') '------------------------'
      do 80 k=1,kmax
        kday = f(k)
        write(3,'(i10,f12.5)') kday,rm(k)
  80  continue
c
      stop
      end
c------------------------------------------------------------------------------
      subroutine meanv(y,n,ymean,kk)
c
c--find and possibly remove mean value from a data set
c--kk = 0 find mean
c--kk = 1 find and remove mean
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      dimension y(1)
c
      ymean = 0.d0
      do 10 i = 1,n
        ymean = ymean + y(i)
  10  continue
      ymean = ymean/float(n)
      if(kk.ne.1) return
c
      do 20 i = 1,n
        y(i) = y(i)-ymean
  20  continue
c
      return
      end
c----------------------------------------------------------------------------
      subroutine fnday(fname,p,fday,ext)
c
c--creates a filename fname from prefix 'p' plus
c--root 'fday' plus extension 'ext'
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      character*40 fname,fnamec
      character*5 fday,fdayc
      character*3 ext,extc
      character*1 p,dum(40),dum3(3),dum5(5)
      equivalence (fnamec,dum),(extc,dum3),(fdayc,dum5)

c--length of prefix is 1, fday is 5 and extension 3, total 10

      fdayc = fday
      extc = ext
      dum(1) = p
      do 10 i = 2,6
        dum(i) = dum5(i-1)
 10   continue
      dum(7) = '.'
      dum(8) = dum3(1)
      dum(9) = dum3(2)
      dum(10) = dum3(3)
      do 20 i = 11,40
        dum(i) = ' '
 20   continue
      fname = fnamec

      return
      end
c---------------------------------------------------------------------
      subroutine fnext2(fname2,fday,ext)
c
c changes the name of file fname2
c gives extension '.ext', where 'ext' is provided on input
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      character*40 fname2,fnamec
      character*5 fday,fdayc
      character*3 ext,extc
      character*1 dum(40),dum3(3),dum5(5)
      equivalence (fnamec,dum),(extc,dum3),(fdayc,dum5)
c
c original length of filename is 'len0'
c find length of extension 'len3' by searching last 4 characters
c
      fnamec = fname2
      len0 = lstr(fnamec,-40)
      len3 = 0
      ltmp = 0
      do 5 i = 1,4
        l = len0-i+1
        if(l.eq.0) goto 6
        if(dum(l).eq.'.') then
          len3 = ltmp+1
          goto 6
        endif
        ltmp = ltmp+1
   5  continue
c
c set fday to first segment name minus the 'g'
c
   6  do 10 i = 1,5
        dum5(i) = dum(i+1)
  10  continue
      fday = fdayc
c
c  add extension '.ext'
c
  20  extc = ext
      len2 = len0-len3
      dum(len2+4) = dum3(3)
      dum(len2+3) = dum3(2)
      dum(len2+2) = dum3(1)
      dum(len2+1) = '.'
      fname2 = fnamec
c
      return
      end
c------------------------------------------------------------------
      subroutine intch(iday,fday)
c
c--convert integer 'iday' to character string 'fday'
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      character*5 fday,fdayc
      character*1 dum(5)
      equivalence (fdayc,dum)
c
      n = iday
      do 10 i = 1,5
        k = n-10*(n/10)
        dum(6-i) = char(k+48)
        n = n/10
 10   continue
      fday = fdayc

      return
      end
c-----------------------------------------------------------------------------
      subroutine chint(fday,kday)
c
c--converts a 5-character string to an integer
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      character*5 fday,fdayc
      character*1 dum(5)
      equivalence (fdayc,dum)
c
      fdayc = fday
      kday = 0
      do 10 i = 1,5
        k = ichar(dum(6-i))-48
        kday = kday+ k*(10**(i-1))
 10   continue
c
      return
      end
c----------------------------------------------------------------------
      subroutine sort(d,t,n)
c
c Subroutine to sort vector t into ascending order of vector d
c Algorithm is based on Shell-Metzner sort (Miller, A.R., 1982. FORTRAN
c programs for scientists and engineers. Sybex)
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      dimension d(1),t(1)
c
      jump = n
  10  jump = jump/2
      if(jump.eq.0) return
      j2 = n-jump
      do 30 j = 1,j2
        i = j
  20    j3 = i+jump
        if(d(i).le.d(j3)) goto 30
        call swap(d(i),d(j3))
        call swap(t(i),t(j3))
        i = i-jump
        if(i.gt.0) goto 20
  30  continue
      goto 10
c
      end
c....................................................
      subroutine swap(a,b)
c
c swaps any two real variables
c
      implicit real*8(a-h,o-y)
c
      hold = a
      a = b
      b = hold
c
      return
      end
c------------------------------------------------------------------------------
      subroutine pfit(x,y,i1,i2,m,a,b)
c
c--fits a polynomial to data in arrays (x(i), y(i), i=i1,i2)
c--the number of coefficients of the polynomial is m (degree = m-1)
c--the arrays must be dimensioned as follows :
c--       x, y : any length
c--       a    : (10,10)
c--       b    : (10)
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      dimension x(1),y(1),
     &          a(10,10),alud(10,10),b(10)
      data np/10/
c
c--zero arrays
c
      do 5 j = 1,10
        b(j) = 0.d0
        do 5 i = 1,10
          a(i,j) = 0.d0
          alud(i,j) = 0.d0
   5  continue
c
c--set up arrays
c
      do 40 j = 1,m
        sumb = 0.d0
        j1 = j-1
        do 10 i = i1,i2
          sumb = sumb+y(i)*(x(i)**j1)
  10    continue
        b(j) = sumb
        do 30 k =1,m
          jk = j+k-2
          suma = 0.d0
          do 20 i = i1,i2
            suma = suma + x(i)**jk
  20      continue
          a(j,k) = suma
  30    continue
  40  continue
c--solve
      call deqn(.false.,a,alud,b,m,np,6,*999)
c
c--coefficients are returned to main program in array b
        return
 999    write(6,'(a)') ' error returned from deqn '
        stop
        end
c-----------------------------------------------------------------------------
      subroutine deqn(zdet,a,alud,b,n,np,iout,*)
c
c  organises equation solution and determinant calculation
c  equations are   a . x = b   on exit solution is in b
c  if zdet=.true. only determinant of matrix a is returned in b(1)
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      parameter(nmax=20)
      real*8 a(np,np),b(np),indx(nmax),alud(np,np),bsave(nmax)
      data big/1.d50/
c
c if zdet=.true. find the determinant of matrix a
c do lu decomposition, note a is replaced
c
      if(zdet) then
        call ludcmp(a,n,np,indx,d,iout,*99)
        do 11 j=1,n
          if(dabs(d).gt.big) then
            write(iout,'(a)') ' error from deqn. determinant too big'
            goto 99
          endif
          d=d*a(j,j)
  11      continue
        b(1)=d
      else
c
c solve linear system, copy a to alud and b to bsave
c
        do 12 j=1,n
          bsave(j)=b(j)
          do 12 i=1,n
            alud(i,j)=a(i,j)
  12        continue
c--lu decomposition returned in alud
        call ludcmp(alud,n,np,indx,d,iout,*99)
c--back substitution destroys b, leaves alud unchanged
        call lubksb(alud,n,np,indx,b)
c--iterative improvement of b requires original a matrix and b vector
        call mprove(a,alud,n,np,indx,bsave,b)
      endif
      return
c
  99  return 1
      end
c---------------------------------------------------------------------------
      subroutine ludcmp(a,n,np,indx,d,iout,*)
c
c lower-upper decomposition
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      parameter (nmax=20)
      real*8 a(np,np),indx(n),vv(nmax)
      data tiny/1.d-20/
c
      d=1.d0
      do 12 i=1,n
        aamax=0.d0
        do 11 j=1,n
          aamax=dmax1(aamax,dabs(a(i,j)))
  11      continue
        if(aamax.eq.0.d0) then
          write(iout,'(a)') ' error from ludcmp. singular matrix '
          return 1
        endif
        vv(i)=1.d0/aamax
  12    continue
      do 19 j=1,n
        if(j.gt.1) then
          do 14 i=1,j-1
            sum=a(i,j)
            if(i.gt.1) then
              do 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
  13            continue
              a(i,j)=sum
            endif
  14      continue
        endif
        aamax=0.d0
        do 16 i=j,n
          sum=a(i,j)
          if(j.gt.1) then
            do 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
   15         continue
            a(i,j)=sum
          endif
          dum=vv(i)*dabs(sum)
          if(dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
  16      continue
        if(j.ne.imax) then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
  17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(j.ne.n) then
          if(a(j,j).eq.0.d0) a(j,j)=tiny
          dum=1.d0/a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
  18        continue
        endif
  19    continue
      if(a(n,n).eq.0.d0) a(n,n)=tiny
c
      return
      end
c----------------------------------------------------------------------------
      subroutine lubksb(a,n,np,indx,b)
c
c backsubstitution
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      real*8 a(np,np),indx(n),b(n)
c
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if(ii.ne.0) then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
  11        continue
        elseif(sum.ne.0.d0) then
          ii=i
        endif
        b(i)=sum
  12    continue
      do 14 i=n,1,-1
        sum=b(i)
        if(i.lt.n) then
          do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
  13        continue
        endif
        b(i)=sum/a(i,i)
  14    continue
c
      return
      end
c-------------------------------------------------------------------
      subroutine mprove(a,alud,n,np,indx,b,x)
c
c improves solution to de's
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      parameter (nmax=20)
      real*8 a(np,np),alud(np,np),indx(n),b(n),x(n),r(nmax)
c
      do 12 i=1,n
        sdp=-b(i)
        do 11 j=1,n
          sdp=sdp+a(i,j)*x(j)
  11      continue
        r(i)=sdp
  12    continue
      call lubksb(alud,n,np,indx,r)
      do 13 i=1,n
        x(i)=x(i)-r(i)
  13    continue
c
      return
      end
c----------------------------------------------------------------------------
      function lstr(name,n)
c
c function to return the length of a string. On input the string
c 'name' has a nominal length n; if n is negative, blank characters
c at the beginning of 'name' are first truncated.
c
      character*1 name(1)
c
c--zero length case first
      if(n.eq.0) then
        lstr = 0
        return
      endif
      k = 0
      ld = n
c
c--ld is negative, we truncate leading blanks
      if(ld.lt.0) then
        ld = -ld
        j = 0
  10    j = j+1
        if(name(j).eq.' ') goto 10
        nb = j-1
c
c--there are nb blanks to be erased in front of 'name'
        if(nb.gt.0) then
          do 20 i = nb+1,ld
            name(i-nb) = name(i)
            name(i) = ' '
  20      continue
c
c--reduce length of string by number of blanks
          ld = ld-nb
        endif
      endif
c
c--ld is non-zero, find length of string
c--first set true length equal to nominal length
      k = ld+1
c--and decrease k until a non-blank character appears or beginning
c--of string is reached
  40  k = k-1
      if(k.gt.0.and.name(k).eq.' ') goto 40
      lstr = k
c
      return
      end