#
# How to use this:
#
# include 'oneoverfnoise.g'
# one := oneoverfnoise();
# one.setup(n=5000, dt=0.01, seed=1005, power=1)
#  #  one.makenoiseguts();
#  #  one.scalenoise(targetmean=1.0, targetrms=0.0001, timerms=1);
# one.makenoise(targetmean=1.0, targetrms=.001, timerms=1);
# myseries := one.getseries();
# gave := one.average(t=5, dt=0.1);
# one.plotseries();
# one.done();
#
include 'randomnumbers.g'
include 'pgplotter.g'
include 'fitting.g'                  
include 'fftserver.g'

oneoverfnoise := subsequence() {

  private := [=];
  private.a := F;   # initial vector of random numbers
  private.freq := F;   # vector of FREQs
  private.timeseries := F;   # final output time series
  private.nx := F;
  private.power := 1.0;
  private.pg := F;
  private.dt := 0.01;   # delta time between gain samples

  self.type := function() {
    return 'oneoverfnoise';
  }

  self.done := function() {
    wider private;
    wider self;

    private.a := F;
    private.freq := F;  
    private.timeseries := F; 
    if (!is_boolean(private.pg)) {
       private.pg.done();
       private.pg := F;
    }
    private := F;
    self := F;
    return T;
  }

# rms (about the average) the time series starting at time = t, going over dt
  self.rms := function(t=0, dt=0.1) {
     if (is_boolean(private.timeseries)) {
	print 'Gotta run makerandom() and makenoise() first!';
	return F;
     }
     nave := max(as_integer(dt / private.dt), 2);

     istart := as_integer( t / private.dt ); # can be 0
     iend := istart + nave - 1;
#     print 'istart = ', istart, '   iend = ', iend,  '   nave = ', nave;
     ave := sum(private.timeseries[([istart:iend] % 
               private.nx) +1] ) / nave;
     rms := sqrt( sum( (private.timeseries[([istart:iend] % private.nx)+1]-ave)^2 )/nave);
     return rms;
  }


# average the time series starting at time = t, going over dt
  self.average := function(t=0, dt=0.1) {
     if (is_boolean(private.timeseries)) {
	print 'Gotta run makerandom() and makenoise() first!';
	return F;
     }
     nave := max(as_integer(dt / private.dt), 1);

     istart := as_integer( t / private.dt ); # can be 0
     iend := istart + nave - 1;
     if (istart < 0 ) {
	istart := istart + private.nx;
	iend   := iend   + private.nx;
     }
#    print 'istart = ', istart, '   iend = ', iend;
     ave := sum(private.timeseries[([istart:iend] % 
                private.nx) +1] ) / nave;
     return ave;
  }

  self.setup := function(n=5000, dt=0.01, seed=1005, power=1.0) {
    wider private;
    rand := randomnumbers();
    private.dt := dt;
    private.seed := seed;
    rand.reseed(private.seed);
    private.nx := n;
    private.power := power;
    private.a := rand.normal(0, 1.0, private.nx);
    rand.done();
    return T;
  }

  self.getseries := function() {
    myvec := private.timeseries;
    return myvec;
  }

  # timerms is the time [s] over which the targetrms is calculated
  self.makenoise := function(targetmean=0.0, 
				targetrms=.0001, 
				timerms=1) {
    wider private;
    self.makenoiseguts();
    self.scalenoise(targetmean, targetrms, timerms);
    return T;
  } 

  # timerms is the time [s] over which the targetrms is calculated
  # WE NEED TO AVERAGE TO THAT TIME!
  self.scalenoise := function(targetmean=1.0, targetrms=0.0001, timerms=1) {
    wider private;

    if (is_boolean(private.timeseries)) {
      print 'Gotta run makerandom() first, man!'
      return F;
    }

    myave := sum(private.timeseries)/len(private.timeseries);
    private.timeseries := private.timeseries - myave;

#  OLD code:
#  print 'myave = ', myave;
#    nsegments := as_integer(private.dt * private.nx / timerms);
#  print 'nsegments = ', nsegments;
#    myrms := [1:nsegments];
#    for (i in [1:nsegments]) {
#       myrms[i] := self.rms(t=i, dt=timerms);
#    }
#    therms := sum( myrms )/nsegments;	
#
#  NEW CODE:

   nn := as_integer(timerms / private.dt);
   newtimeseries := private.averesample( private.timeseries, nn );
   therms := private.stepvar(newtimeseries, 1);
#  print 'therms = ', therms;
    private.timeseries := private.timeseries * targetrms / therms + targetmean;
    return T;
  }

# calculate the allan standard deviation
#
   self.allansd := function(ref y=[], dt=0.1) {

    xalvar1 := as_integer(1.3 ^[1:30]);
    xalvar := unique(xalvar1[ (xalvar1 < len(y)/5 ) ]);
    yalvar := 0*xalvar;

    for (i in ind(xalvar)) {
     newvec := private.averesample(y, xalvar[i]);
     yalvar[i] := private.stepvar(newvec) ;
     print (dt * xalvar[i]), yalvar[i]
    }

   results := [=];
   results.x := xalvar;
   results.y := yalvar;
   return results;
  } 

  self.makenoiseguts := function() {
    wider private;

    if (is_boolean(private.a)) {
      print 'Gotta run makerandom() first, man!'
      return F;
    }
    server := fftserver()
    b := server.realtocomplexfft(private.a);
    private.freq :=  private.nx/2 - [0:(private.nx-1)]
    private.freq[private.nx/2+1] := 1.0;
    oneoverf := 1.0/abs(private.freq);
    oneoverf2 := (oneoverf)^(private.power/2)
    oneoverf2[private.nx/2+1] := 0.0;
    b := b * oneoverf2;
    server.complexfft(b, 1);
    server.done()
    private.timeseries := real(b);
    b := F;
    return T;
  }

  self.plotseries := function() {
    wider private;
    wider self;

    if (is_boolean(private.timeseries)) {
	print 'Gotta run makerandom() and makenoise() first!';
	return F;
    }
    if (is_boolean(private.pg)) {
       private.pg := pgplotter();
    }
    private.pg.page();
    private.pg.env(min(private.freq), max(private.freq), 
	min(private.timeseries), max(private.timeseries), 0, 0);
    private.pg.line(private.freq, private.timeseries);
    return T;
  }


#  private functions:

#  resamples and averages
private.averesample := function( ref vec=[], nsamp ) {
   nnewvec := as_integer( len(vec)/nsamp );
   newvec := 0 * [1:nnewvec];   
   for (i in [1:nsamp]) {
     iii := ind(newvec)
     newvec [ iii ] +:= vec[ (nsamp * (iii-1) + i ) ]
   }
   newvec /:= nsamp;
   return ref newvec
}

private.stepvar := function( ref vec=[], step=1) {
   nn := len(vec);
   if (nn <= (step+2)) {
         print 'Not enough elements in stepvar()';
	 return F;
   }
   var := sqrt( sum( (vec[1:(nn-step)] - vec[(1+step):nn])^2 )/(nn-step-1) );
   return var;
}


  
}

# end of product definition

#
#  make a log log plot of the power spectrum of timeseries
#
#server := fftserver()
#b := server.realtocomplexfft(timeseries);
#b := (abs(b))^2
#server.done()
#pg.done();
#pg := pgplotter()
#pg.env(0.0, 2.7, 0, 10, 0, 30)
#lgfreq := log( freq[1:nxo2] )
#lgb    := log( b[1:nxo2] )
#pg.pnts(lgfreq, lgb, 2)
#
#

testallansd := function(seed=1005, powerlaw=1.0, amp=0.001) {

 one := oneoverfnoise();
 one.setup(n=1000, dt=0.1, seed=seed, power=powerlaw)
 one.makenoise(targetmean=1.0, targetrms=amp, timerms=1);
 myseries := one.getseries();
 rec := one.allansd( myseries, dt=0.1 );
 one.done();
 return rec;
}



testoneoverf := function(seed=1005, power=1) {

 include 'oneoverfnoise.g'
 one := oneoverfnoise();
 one.setup(n=5000, dt=0.01, seed=seed, power=power)
  #  one.makenoiseguts();
  #  one.scalenoise(targetmean=1.0, targetrms=0.0001, timerms=1);
 one.makenoise(targetmean=1.0, targetrms=.0001, timerms=1);
 # noisevec := one.getseries( );
 gave := one.average(t=5, dt=0.1);
 print 'Average gain = ', gave;
 one.plotseries();
 one.setup(n=5000, dt=0.01, seed=1005, power=power)
 one.done();
 return T;

}


testaverage := function(seed=1005, power=1) {

 seed := 1005;  power := 1;
 include 'oneoverfnoise.g'
 one := oneoverfnoise();
 n:= 2^14;
 one.setup(n=n, dt=0.01, seed=seed, power=power)
 one.makenoise(targetmean=1.0, targetrms=.0001, timerms=1);
 gainvec := one.getseries();

 myrand := randomnumbers();
 myrand.reseed( seed +100 );
 thermalvec := myrand.normal(0, 1.0, n);
 myrand.done();
 one.done();

 avedown := [1:14];
 ndown   := [1:14];
 for (i in ndown) {
   gain := gainvec[1:(2^i)];
   thermal := thermalvec[1:(2^i)];
   output := gain * thermal;
   avedown[i] := sqrt( sum(output^2 ) )/ len(output)
 }

 logavedown := log(avedown)
 pg := pgplotter();
    pg.page();
    pg.env(min(ndown), max(ndown), 
	min(logavedown), max(logavedown), 0, 20);
    pg.pnts(ndown, logavedown, 2);

    myfit:=fitter();                  
    myfit.init(n=2);                  
    myfit.makepoly(x,y);              
    myfit.fit();                      
    sol:=myfit.solution();            
    myfit.done();

 return T;

}