program noisecon
      dimension per(21),a(21),b(21)
      real noise,noisin,noisep,noised,lnm,lnmd,mil
      integer ex,dva
      character dis,amp,rms,ban,okt,fre,whatever*24,ok,amp1*13
c  
c   USGS New Low Noise Model (Peterson 1993)
c  
      data per(1),a(1),b(1)    / 0.1,-162.36,5.64 /
      data per(2),a(2),b(2)    / 0.17,-166.7,0 /
      data per(3),a(3),b(3)    / 0.4,-170,-8.3 /
      data per(4),a(4),b(4)    / 0.8,-166.4,28.9 /
      data per(5),a(5),b(5)    / 1.24,-168.6,52.48 /
      data per(6),a(6),b(6)    / 2.4,-159.98,29.81 /
      data per(7),a(7),b(7)    / 4.3,-141.1,0 /
      data per(8),a(8),b(8)    / 5,-71.36,-99.77 /
      data per(9),a(9),b(9)    / 6,-97.26,-66.49 /
      data per(10),a(10),b(10) / 10,-132.18,-31.57 /
      data per(11),a(11),b(11) / 12,-205.27,36.16 /
      data per(12),a(12),b(12) / 15.6,-37.65,-104.33 /
      data per(13),a(13),b(13) / 21.9,-114.37,-47.1 /
      data per(14),a(14),b(14) / 31.6,-160.58,-16.28 /
      data per(15),a(15),b(15) / 45,-187.5,0 /
      data per(16),a(16),b(16) / 70,-216.47,15.7 /
      data per(17),a(17),b(17) / 101,-185,0 /
      data per(18),a(18),b(18) / 154,-168.34,-7.61 /
      data per(19),a(19),b(19) / 328,-217.43,11.9 /
      data per(20),a(20),b(20) / 600,-258.28,26.6 /
      data per(21),a(21),b(21) / 10000,-346.88,48.75 /
c
      pi=4.*atan(1.)
      write(6,*)
      write(6,*) ' *** CONVERSION OF NOISE SPECTRA ***',
     &'    Adapted from the GFA-Basic version: Apr 1999, EW'
      write(6,*)
c
c   Definition of input units
c     write(6,*)
      write(6,'(a/a)') ' Please describe your input data by answering',
     &' multiple-choice questions with initial of keyword'
      write(6,'(/a/)') ' Fragen mit dem Anfangsbuchstaben des Stichworts
     & beantworten!'

c questions:
    1 write(6,*)
    2 write(6,*) ' spectrum of Displacement, Velocity or Acceleration ?'
      read(5,*) dis
      call up(dis)
      if(.not.(dis.eq.'D'.or.dis.eq.'V'.or.dis.eq.'A')) goto 2
    3 write(6,*) ' scaled as Amplitude, Power, Decibels, or Millimeters
     &from baseline of plot ?'
      read(5,*) amp
      call up(amp)
      if(.not.(amp.eq.'A'.or.amp.eq.'P'.or.amp.eq.'D'.or.amp.eq.'M')) 
     &  goto 3
      IF (amp.eq.'A') THEN
    4 write(6,*) ' Rms, Average peak, Zero-to-peak or Peak-to-peak ?'
      read(5,*) rms
      call up(rms)
      if(.not.(rms.eq.'R'.or.rms.eq.'A'.or.rms.eq.'Z'.or.rms.eq.'P'))
     &  goto 4
      ELSE IF (amp.eq.'M') then
      write(6,*) ' how many millimeters per 20 dB ?'
      read(5,*) mil
      ENDIF
    5 IF (amp.eq.'A') THEN
      write(6,*) ' length unit is 10**-n meters, enter n: '
      read(5,*) ex
      ELSE IF (amp.eq.'P') then
      write(6,*) ' square of length unit is (10**-n meters)**2, enter n:
     & '
      read(5,*) ex
      ELSE IF (amp.eq.'D') then
      write(6,*) ' reference is [square of] 10**-n meters, enter n: '
      read(5,*) ex
      ELSE IF( amp.eq.'M') then
      write(6,*) ' level of baseline is [square of] 10**-n meters, enter 
     &n:'
      read(5,*) ex
      ENDIF
      if(.not.(ex.ge.0.and.ex.le.12)) goto 5
    6 write(6,*) ' Constant or Proportional bandwidth ?'
      read(5,*) ban
      call up(ban)
      if(.not.( ban.eq.'C'.or.ban.eq.'P')) goto 6
      IF (ban.eq.'C') THEN
      write(6,*) ' bandwidth in Hz ?'
      read(5,*) bw
      ELSE
    7 write(6,*) ' is the bandwidth specified as fraction of Octave or D
     &ecade ?'
      read(5,*) okt
      call up(okt)
      if(.not.(okt.eq.'O'.or.okt.eq.'D')) goto 7
      write(6,*) ' what fraction (enter integer denominator)'
      read(5,*) bw
      ENDIF
    8 write(6,*) ' will you enter Frequency (Hz) or Period (s) ?'
      read(5,*) fre
      call up(fre)
      if(.not.(fre.eq.'F'.or.fre.eq.'H'.or.fre.eq.'P'.or.fre.eq.'S')) 
     &  goto 8
      write(6,*)
      write(6,*)
      write(6,*) ' ===>  You say that our input data are in units of:'
      write(6,*)
      IF (amp.eq.'D') THEN
        write(6,*) ' decibels relative to'
      ELSE IF (amp.eq.'M') then
        write(6,*) 20/mil,' decibels relative to'
      ENDIF
      ibw=nint(bw)
      IF (ban.eq.'C') THEN
        write(whatever,*) bw,' Hz'
      ELSE IF (okt.eq.'O') then
        write(whatever,*) '1/',ibw,' octave'
      ELSE IF (okt.eq.'D') then
        write(whatever,*) '1/',ibw,' decade'
      ENDIF
      IF (dis.eq.'D') THEN
        dva=0
      ELSE IF (dis.eq.'V') then
        dva=1
      ELSE IF (dis.eq.'A') then
        dva=2
      ENDIF
      IF (amp.eq.'A') THEN
        IF (rms.eq.'R') THEN
          amp1=' rms         '
        ELSE IF (rms.eq.'A') then
          amp1=' average peak'
        ELSE IF (rms.eq.'Z') then
          amp1=' zero-to-peak'
        ELSE IF (rms.eq.'P') then
          amp1=' peak-to-peak'
      ENDIF
      write(6,*) ' (10**-',ex,' m/s**',dva,amp1,') per square root of'
     &, whatever
      ELSE
        write(6,*) ' (10**-',ex,' m/s**',dva,' ) squared per',whatever
      ENDIF
      write(6,*)
      write(6,*) ' OK - Yes or No ?'
      read(5,*) ok
      call up(ok)
      IF(ok.ne.'Y') goto 1
c
c   Factor for relative bandwidth
c
      bw6=10**(1./12.)-10**(-1./12.)
      IF (ban.eq.'P') THEN
        IF (okt.eq.'O') THEN
          basis=2
        ELSE
          basis=10
        ENDIF
        bw=basis**(0.5/bw)
        bw=bw-1/bw
      ENDIF
      icount=0

c  
c   Conversion of input data
c  
   10 icount=icount+1
      IF (mod(icount,5).eq.1) THEN
        write(6,*)
        write(6,*) ' Noise power densities in dB re (1 m/s**2)**2/Hz.', 
     &  ' Noise amplitudes in dB re.'
        write(6,*) ' 1 m/s**2 rms in 1/6 decade, equiv. to avg. pk.',
     &  ' amplitude in 1/3 octave'
        write(6,*) ' LNM (New Low Noise Model) is given in your units,'
     &,' power density, and amplitude'
      ENDIF
      write(6,*)
   11 IF (fre.eq.'F'.or.fre.eq.'H') THEN
        write(6,*) ' Enter frequency, noise (0,0=end)'
        read(5,*) f,noise
        if(f.eq.0.) stop
        p=1/f
      ELSE
        write(6,*) ' Enter period, noise (0,0=end)'
        read(5,*) p, noise
        if(p.eq.0.) stop
        f=1/p
      ENDIF
      if((amp.eq.'A'.or.amp.eq.'P').and.noise.eq.0.) then
        write(6,*) 'illegal input for noise'
        goto 11
      endif
      noisin=noise
      omega=2*PI*f
c   absolute bandwidth of input data
      IF (ban.eq.'P') THEN
        bwin=bw*f
      ELSE
        bwin=bw
      ENDIF
      bwout=bw6*f
c   Low Noise Model for comparison
      lnmd=fnlnm(p,per,a,b)
      lnm=lnmd+10*LOG10(bwout)
c   conversion from amplitude to power
      IF (amp.eq.'A') THEN
        IF (rms.eq.'A') THEN
          noise=noise/1.25
        ELSE IF (rms.eq.'Z') then
          noise=noise/3
        ELSE IF (rms.eq.'P') then
          noise=noise/6
        ELSE IF (rms.ne.'R') then
          write(6,*) ' Illegal specification of units'
          STOP
        ENDIF
        noisep=noise**2
c   conversion from decibels to power
      ELSE IF (amp.eq.'D') then
        noisep=10**(noise/10)
      ELSE IF (amp.eq.'M') then
        noisep=10**(noise/mil*2)
      ELSE IF (amp.eq.'P') then
        noisep=noise
      ELSE
        write(6,*) ' Illegal specification of units'
        STOP
      ENDIF
c   conversion to acceleration
      noisep=noisep*omega**(4-2*dva)
c   acceleration power density
      noised=noisep/bwin
c   conversion to decibels of rms amplitude and power density
      ampl=10*LOG10(noised*bwout)-20*ex
      dens=10*LOG10(noised)-20*ex
c   Convert Low Noise Model into input units
      IF (amp.eq.'A') THEN
        noisin=noisin*10**((lnmd-dens)/20)
      ELSE IF (amp.eq.'D') then
        noisin=noisin+lnmd-dens
      ELSE IF (amp.eq.'M') then
        noisin=noisin+(lnmd-dens)/mil
      ELSE IF (amp.eq.'P') then
        noisin=noisin*10**((lnmd-dens)/10)
      ELSE
        write(6,*) ' Illegal specification of units'
        STOP
      ENDIF
      write(6,*) ' frequency    period       LNM    LNM pd   your pd   L
     &NM amp  your amp  diff'
      write(6,101) f,p,noisin,lnmd,dens,lnm,ampl,dens-lnmd
  101 format(3f10.3,4f10.0,f6.0)
      goto 10
      END

c   evaluation of low noise model
      FUNCTION fnlnm(p,per,a,b)
      dimension per(21),a(21),b(21)
      IF (p.lt.0.1) THEN
        write(6,*) ' NLNM undefined, OLNM used'
        fnlnm=-168.
      ELSE IF (p.le.100000.) then
        do k=1,20
          IF (p.lt.per(k+1)) goto 20
        enddo
   20   fnlnm=a(k)+b(k)*LOG10(p)
      ELSE
        write(6,*) ' NLNM undefined'
        fnlnm=0.
      ENDIF
      END

      subroutine up(ch)
      character ch
      if(iachar(ch).gt.96) ch=achar(iachar(ch)-32)
      return
      end