#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>

float nlnm();

main (arc,argv)
int     arc;
char    *argv[];
{

   int True, dva, icount, ibw;
   div_t mod;
   char STRING[129], whatever[80], amp1[80];
   char dis, amp, rms, ban, oct, fre;
   float mil, ex, bw, bw6, basis, f, p, noise, noisin, omega, bwin, bwout;
   float PI, lnm, lnmd, noisep, noised, ampl, dens, unsinn;
   
   static float per[21]= { 0.1,0.17,0.4, 0.8,1.24,2.4,4.3,5,
                         6,10,12,15.6,21.9,31.6,45,70,101,
                        154,328,600,10000
   };
   static float a[21]=   {-162.36,-166.7,-170,-166.4,-168.6,-159.98,-141.1,-71.36,
                         -97.26,-132.18,-205.27,-37.65,-114.37,-160.58,-187.5,-216.47,-185,
                         -168.34,-217.43,-258.28,-346.88
   };
   static float b[21]=   {5.64,0,-8.3,28.9,52.48,29.81,0,-99.77,
                        -66.49,-31.57,36.16,-104.33,-47.1,-16.28,0,15.7,0,
                         -7.61,11.9,26.6,48.75
   };
   PI= atan(1.0)*4.0;
   unsinn=777.0e24;
   printf("\n *** CONVERSION OF NOISE SPECTRA ***    E. Wielandt 1995\n");
   printf("\n *** New Low Noise Model by Jon Peterson, USGS Open-File Report 93-322\n");
   printf("\n *** C version by M. Baer, 1996\n");
   printf("\n *** ");

/****************************
' Definition of input units
****************************/

questions:   printf("please describe your input data by answering multiple-choice\n");
   printf("     questions with initial of keyword\n");
  printf("\n ");

   True=1;
   while ( True ) {
      printf("spectrum of Displacement, Velocity or Acceleration [D/V/A] \n");
      fflush(stdin);
      fgets(STRING,128,stdin);
      STRING[0]= toupper(STRING[0]);
      switch (STRING[0])
      {
	 case 'D':
	 case 'V':
	 case 'A':
	    dis=STRING[0];
	    True= 0;
	    break;
      }
   }
   True=1;
   while ( True ) {
      printf("scaled as Amplitude, Power, Decibels,\n          Millimeters from baseline of plotor [A/P/D/M] \n");
      fflush(stdin);
      fgets(STRING,128,stdin);
      STRING[0]= toupper(STRING[0]);
      switch (STRING[0])
      {
	 case 'D':
	 case 'P':
	 case 'A':
	 case 'M':
	    amp=STRING[0];
	    True= 0;
	    break;
      }
   }
   if ( amp == 'A' ) {
      True=1;
      while ( True ) {
         printf("Rms, Average peak, Zero-to-peak or Peak-to-peak [R/A/Z/P] \n");
         fflush(stdin);
         fgets(STRING,128,stdin);
         STRING[0]= toupper(STRING[0]);
         switch (STRING[0])
         {
	    case 'R':
	    case 'P':
	    case 'A':
	    case 'Z':
	       rms=STRING[0];
	       True= 0;
	       break;
         }
      }
   }
   if ( amp == 'M' ) {
      printf("how many millimeters per 20 dB\n");
      fflush(stdin);
      fscanf(stdin,"%f",&mil);
   }
   ex= -1;
   while ( ex < 0 || ex > 12 ) {
      switch (amp)
      {
	 case 'A':
	    printf("length unit is 10^-n meters, enter n: \n");
	    break;
	 case 'P':
	    printf("square of length unit is (10^-n meters)^2, enter n: \n");
	    break;
	 case 'D':
	    printf("reference is [square of] 10^-n meters, enter n: \n");
	    break;
	 case 'M':
	    printf("level of baseline is [square of] 10^-n meters, enter n: \n");
	    break;
      }
      fflush(stdin);
      fscanf(stdin,"%f",&ex);
   }
   True=1;
   while ( True ) {
      printf("Constant or Proportional bandwidth [C/P] \n");
      fflush(stdin);
      fgets(STRING,128,stdin);
      STRING[0]= toupper(STRING[0]);
      switch (STRING[0])
      {
	 case 'C':
	 case 'P':
	    ban=STRING[0];
	    True= 0;
	    break;
      }
   }
   if ( ban == 'C' ) {
      printf("bandwidth in Hz\n");
      fflush(stdin);
      fscanf(stdin,"%f",&bw);
   }
   else {
      True=1;
      while ( True ) {
         printf("is the bandwidth specified as fraction of Octave or Decade [O/D] \n");
         fflush(stdin);
         fgets(STRING,128,stdin);
         STRING[0]= toupper(STRING[0]);
         switch (STRING[0])
         {
	    case 'O':
	    case 'D':
	       oct=STRING[0];
	       True= 0;
	       break;
         }
      }
      printf("what fraction (enter integer denominator)\n");
      fflush(stdin);
      fscanf(stdin,"%f",&bw);      
   }
   True=1;
   while ( True ) {
      printf("will you specify Frequency (Hz) or Period (s) [F/H/P/S] \n");
      fflush(stdin);
      fgets(STRING,128,stdin);
      STRING[0]= toupper(STRING[0]);
      switch (STRING[0])
      {
	 case 'F':
	 case 'H':
	 case 'P':
	 case 'S':
	    fre=STRING[0];
	    True= 0;
	    break;
      }
   }
   printf("\n\n===>  You say that your input data are in units of: \n\n");
   if ( amp == 'D' ) {
      printf("decibels relative to:\n");
   }
   else if ( amp == 'M' ) {
      printf("20/%f decibels relative to:\n",mil);
   }
   switch (ban) {
      case 'C':
	 sprintf(whatever,"%f Hz ",bw);
	 break;
      case 'P':
      ibw=bw+0.001;
      switch (oct) {
      case 'O':
	 sprintf(whatever,"1/%d octave ",ibw);
	 break;
      case 'D':
	 sprintf(whatever,"1/%d decade ",ibw);
	 break;
   }
         break;
   }
   switch (dis) {
      case 'D':
	 dva= 0;
	 break;
      case 'V':
	 dva= 1;
	 break;
      case 'A':
	 dva= 2;
	 break;
   }
   if ( amp == 'A') {
      switch (rms) {
         case 'R':
	    strcpy(amp1," rms ");
	    break;
         case 'A':
	    strcpy(amp1," average peak ");
	    break;
         case 'Z':
	    strcpy(amp1," zero-to-peak ");
	    break;
         case 'P':
	    strcpy(amp1," peak-to-peak ");
	    break;
      }
      printf("( 10^-%d m/s^%d%s) per square root of %s\n",(int)ex,dva,amp1,whatever);
   }
   else {
      printf("( 10^-%d m/s^%d) squared per %s\n",(int)ex,dva,whatever);
   }
   printf("\nOK - yes or no [Y/N]:\n");
   fflush(stdin);
   fgets(STRING,128,stdin);
   STRING[0]= toupper(STRING[0]);
   if ( STRING[0] == 'N' ) goto questions;
   /******* Factor for relative bandwidth *********/
   bw6=pow(10.0,1.0/12.0) - pow(10.0,-1.0/12.0);
   /*bw6=10^(1/12)-10^(-1/12)*/
   if ( ban == 'P' ) {
      if ( oct == 'O' ) {
         basis= 2;
      }
      else {
	 basis= 10;
      }
      bw=pow(basis,0.5/bw);
      bw-= 1./bw;
   }
  icount=0;
  True=1;
  while (True)
  {
      icount++;
      mod= div(icount,5);
      if (mod.rem == 1 ) {
	 printf("\n");
         printf(" Noise power densities in dB re (1 m/s^2)^2/Hz ; Noise amplitudes in dB re.\n");
         printf(" 1 m/s^2 rms in 1/6 decade, equiv. to avg. pk. amplitude in 1/3 octave\n");
         printf(" LNM (New Low Noise Model) is given in your units, power density, and amplitude\n");
      }
      printf("\n");
      noise=unsinn;
      if ( fre == 'F' || fre == 'H' )
 { printf("frequency, noise level (separated by <space>; 0 0 = end):");
         fflush(stdin);
	 fscanf(stdin,"%f %f",&f,&noise);
	 if ( f == 0.0 ) exit(0);
	 p=1.0/f;
      }
      else {
	 printf("period, noise level ( 0 0 = end):");
         fflush(stdin);
	 fscanf(stdin,"%f %f",&p,&noise);
	 if ( p == 0.0 ) exit(0);
	 f=1.0/p;
      }
      if (noise == unsinn) {
        printf("!!! use <space> as a separator !!!\7\n");
      }
      noisin=noise;
      omega=2*PI*f;
      /** absolute bandwidth of input data **/
      if ( ban == 'P' ) {
	 bwin=bw*f;
      }
      else {
	 bwin=bw;
      }
      bwout=bw6*f;
      /** Low Noise Model for comparison **/
      lnmd= nlnm(p,per,a,b);
      lnm= lnmd+10.0*log10(bwout);
      /** conversion from amplitude to power **/
      switch (amp) {
	    case 'A': {
	       switch (rms) {
	          case 'A':
	             noise/=1.25;
	             break;
	          case 'Z':
	             noise/=3;
	             break;
	          case 'P':
	             noise/=6;
	             break;
                  case 'R':
                     break;
	          default:
		     printf("Illegal specification of units\n");
		     exit;
	       }
	       noisep= noise*noise;
	       break;
	    }
	    case 'D':
	       noisep=pow(10.0,noise/10.0);
	       break;
	    case 'M':
	       noisep=pow(10.0,noise*2.0/mil);
	       break;
	    case 'P':
	       noisep=noise;
	       break;
	    default: 
	       printf("Illegal specification of units\n");
	       exit;
      }
      /** conversion to acceleration **/
      noisep=noisep*pow(omega,(float)(4-2*dva));
      /**  acceleration power density **/
      noised=noisep/bwin;
      /**  conversion to decibels of rms amplitude and power density **/
      ampl=10.0*log10(noised*bwout)-20.0*ex;
      dens=10.0*log10(noised)-20.0*ex;
      /**  Convert Low Noise Model into input units **/
      switch (amp) {
	 case 'A':
	    noisin*= pow(10.0,(lnmd-dens)/20);
	    break;
	 case 'D':
	    noisin+= lnmd-dens;
	    break;
	 case 'M':
	    noisin+= (lnmd-dens)/mil;
	    break;
	 case 'P':
	    noisin*=pow(10.0,(lnmd-dens)/10);
	    break;
	 default: 
	    printf("Illegal specification of units\n");
	    exit;
      }
      printf("frequency   period      LNM   LNM pd  your pd   LNM amp  your amp  diff\n");
      printf("%8.3f %8.3f %8.3f%9.0f %8.0f %8.1f %8.1f %8.1f\n",f,p,noisin,lnmd,dens,lnm,ampl,dens-lnmd);
  }
}
/********************************************************/
/********
' evaluation of low noise model
********/
float nlnm(p,per,a,b)
   float p;
   float *per,*a,*b;
{
   int i,k;
   if ( p < 0.1 ) {
      printf(" No NLNM available, OLNM used\n");
      return(-168.0);
   }
   if ( p > 1000000 ) {
      printf(" No NLNM available\n");
      return (0.0 );
   }
   for (k=1,i=0;k<21;k++,i++) {
      if ( p < per[k] ) break;     
   }
   return (a[i]+b[i]*log10(p));
}