# phasemonitor.g   incapsulates the behavior of the
#                  300m baseline site testing interferometer


# pragma include once;

include 'fitting.g'
include 'image.g'

phasemonitordemo := function(atmos='ATMOS', oversample=15) {


  im := image(atmos);

  phmon	    := phasemonitor();
  phmon.setsite('alma');

  cs := im.coordsys();
  imdelta := cs.i()[1];
  cs.done();
  arr := im.getchunk();
  im.done();

  phmon.setarray( arr, imdelta );
  arr := F;
  rec := phmon.getspatialstructurefunc(doplot=T, oversample=oversample)
  print rec;
  phmon.done();
  return rec;
}




phasemonitor := subsequence() {


  private := [=];

  private.arr       := F;               # atmos array
  private.imdelta   := 0.5;             # cellsize
  private.shape     := [];		# shape of array

  private.baseline  := 300;         # site testing interferometer bl [m]
  private.freq      := 11.198e+9;   # sti freq [Hz]
  private.c         := 3e+8;        # [m/s]
  private.factor    := 6.3/1000;    # path [m] caused by 1mm PWV
  private.el        := 36.0    ;    # elevation angle [deg]

  private.pg := F;		    # pg plotter tool


self.type := function() {
   note ('I simulate phase monitor functions');
   return 'phasemonitor';
}

self.getprivate := function() {
        return ref private;
}

self.done := function() {
        wider self, private;
        private.donepgplotter();
        self := F;
        val public := F;
        return T;
}


self.setarray := function(ref arr=[], cell=0.1) {
	wider private;
	private.arr := ref (arr);
	private.shape := shape(arr);
	private.imdelta := cell;
	return T;
}

self.setparms := function(baseline=300, freq=11.198e+9, el=36) {
	wider private;
	  private.baseline  := baseline;
	  private.freq      := freq;
	  private.el        := el;
	  return T;
}
self.setsite := function(site='alma') {
	wider private;

	if (site == 'alma') {
	  private.baseline  := 300;
	  private.freq      := 11.198e+9;
	  private.el        := 36.0;
	  return T;
	} else if (site == 'maunakea') {
	  private.baseline  := 300;
	  private.freq      := 11.7269e+9;
	  private.el	    := 29.0;
	  return T;
	} else if (site == 'vla') {
	  private.baseline  := 300;
	  private.freq      := 11.7269e+9;
	  private.el	    := 49.0;
	  return T;
	} else {
	  print 'Unknown site given: ', site;
	  print 'Known sites are: alma, maunakea, vla';
	  return F;
	}
}




# returns the data for the spatial structure function
# and the fit parameters in a record
self.getspatialstructurefunc := function(doplot=T, oversample=15) {

   results := [=];

   if (is_boolean(oversample)) {
       oversample := 15;
   }
   if (!private.checkarr()) { fail; }

   bmax := as_integer(private.shape[1] / oversample);
   divideby := 1.4
   bmin := as_integer(bmax/divideby);
   nn := 1;
   while (bmin > 2) {
      nn +:= 1;
      bmin /:= divideby;
   }
   results.xx := [1:nn];
   results.rms := results.xx * 0;
   results.xx[nn] := bmax;
   for (i in [(nn-1):1]) {
      results.xx[i] := as_integer(results.xx[i+1]/divideby);
   }

   ny := private.shape[2];
   iy1 := as_integer(0.25 * ny);
   iy2 := as_integer(0.50 * ny);
   iy3 := as_integer(0.75 * ny);

   strip1 := private.arr[,iy1,,,];
   strip2 := private.arr[,iy2,,,];
   strip3 := private.arr[,iy3,,,];
   nstrip := len(strip1);
   for (i in ind(results.xx)) {
     
	one := strip1[(1+results.xx[i]):nstrip];
	other := strip1[1:(nstrip-results.xx[i])];
	diff := one - other;
	ave1 :=  sum( diff^2 )/len(diff) ;

	one := strip2[(1+results.xx[i]):nstrip];
	other := strip2[1:(nstrip-results.xx[i])];
	diff := one - other;
	ave2 :=  sum( diff^2 )/len(diff) ;

	one := strip3[(1+results.xx[i]):nstrip];
	other := strip3[1:(nstrip-results.xx[i])];
	diff := one - other;
	ave3 :=  sum( diff^2 )/len(diff) ;

	results.rms[i] := sqrt( (ave1 + ave2 + ave3)/3.0 );        
   }
   one    := F;
   other  := F;
   diff   := F;
   strip1 := F; strip2 := F; strip3 := F;
   results.xx := results.xx * private.imdelta;

   logxx := log(results.xx);
   logyy := log(results.rms)
   fitresults := [=];
   fitresults := private.lsq (logxx, logyy);
   modyy := fitresults.b + fitresults.a * logxx;
   results.fitamp := 10.0^fitresults.b;
   results.fitexp := fitresults.a;
   fitresults := F;

   if (doplot) {
     private.getpgplotter();
     private.pg.page(); 
     private.pg.sci(1);
     private.pg.env(log(0.8*min(results.xx)), log(1.2*max(results.xx)), 
	log(0.8*min(results.rms)), log(1.2*max(results.rms)), 0, 30)
     private.pg.lab('r [m]', 'rms pwv [mm]', 'PWV Structure Function');
     private.pg.line (logxx, logyy);
     private.pg.sci(2);
     private.pg.line (logxx, modyy);
   }
# calculate the fluctuation level we would have for an interferometer

   dpwv := results.fitamp * private.baseline^results.fitexp;
   dpath := private.factor * dpwv;
   wavelength := private.c / private.freq;
   phase_deg := 360.0 * dpath/wavelength;
   avepwv := sum(private.arr) / len(private.arr);

   results.avepwv := avepwv;
   results.stiphase := phase_deg;
   return results;
} 

########### private functions : members only, please ########

 private.getpgplotter := function() 
  {
    wider self, private;
    if (!is_tool(private.pg)) {
       private.pg := pgplotter();
    }
    return T;
  }

  private.donepgplotter := function() 
  {
    wider self, private;
    if (is_tool(private.pg) ) {
        private.pg.done();
    }
    private.pg := F;
    return T;
  }

  private.checkarr := function()
  {
	if (len(private.arr) <=1 ) {
		print 'phasemonitor::checkarr(): arr not found!';
		return F;
	}
	if (len(private.shape) <= 1) {
		print 'phasemonitor::checkarr(): arr not 2d';
		return F;
	}
	return T;
  }

  private.lsq := function(xx=[], yy=[]) {
    myfit := fitter();
    myfit.init(2);
    myfit.makepoly(xx, yy);
    myfit.fit();
    results := [=];
    results.a :=  myfit.solution()[2];
    results.b :=  myfit.solution()[1];
    myfit.done();
    return results;
  }




}
# end of constructor