# pragma include once;
#
##  2003dec13   LAMA Memo 803 writeup version

include 'quanta.g'
include 'measures.g'
include 'pgplotter.g'
include 'almatau.g'

#const
almasensitivity := subsequence () {

  private := [=];

  private.nants := 64;
  private.diam  := 12.0;	# Diam in meters
  private.baseline:= 12.0;	# average baseline, for resolution
  private.freq := F;		# frequency in GHz

  private.pwv := F; 		# PWV            \____ two words for the same 
  private.tau225 := F;		# Tau at 225 GHz /     concept, in different languages
  private.tatmos := 250;	# atmospheric temperature
  private.tamb   := 270;	# ambient temperature
  private.taufreq := F;
  private.tcmb   := 2.7;	# microwave background
  private.airmass := 1.0;	
  private.almatau := almatau();	# alma tau spectrum
  private.gaincorrect := T 	# increase noise by e^{tau.airmass} ??

  private.etaa := 0.75;		# aperture efficiency
  private.etal := 0.95;		# fprward spillover efficiency ??
  private.etaq := 0.95;		# quantization efficiency

  private.npol := 2;
  private.bandwidth := 8;
  private.freq := 300;
  private.dt := 1;

  private.h := 6.6262e-34;
  private.k := 1.38e-23;
  private.jy := 1e+26;
  private.sensunits := 'Jy';	# J == janskies
# private.sensunits := 'K ';	# K == kelvins

# receiver band information;  note, trx is for SSB, over 80% of the band
# Taken from  May 2001 Specification for Front End Assembly, Wild & Payne
  private.bands    :=   [ 1,   2,   3,  4,   5,   6,   7,   8,   9,   10];
  private.bandlow  :=   [31.3, 67,  84, 125, 163, 211, 275, 385, 602, 787];
  private.bandhigh :=   [45,   90, 116, 163, 211, 275, 370, 500, 720, 950];
  private.bandtrxspec:= [15,   28,  34,  47,  60,  77, 133, 181, 335, 438]; 
#  private.bandtrxspec:= [15,   28,  25.4,  47,  60,  77, 133, 181, 298.4, 438];
  private.bandtrxgoal:= [10,   16,  19,  26,  32,  40,  69,  93, 202, 351]; 
  private.band := F;
  private.rxspec := T;    # else use goal
#
  private.silent := F;
#
# functions
#

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

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

self.silent := function(silent=T) {
  wider private;

  private.silent := silent;

  return T;
}

self.setgaincor := function(doit=T) {
   wider private;
   private.gaincorrect := doit;

   if (! private.silent) {
     if (private.gaincorrect) {
       print 'We will increase the noise by e^{tau.airmass}'
     } else {
       print 'Noise will NOT be increased by  e^{tau.airmass}'
     }
   }
   return T;
}

# automatically calculates aperture efficiency from freq
self.setefficiencies := function(etal=0.95, etaq=0.95) {
  wider private;
  private.etaq := etaq;
  private.etal := etal;
  dummy := self.getetaa();
  return;
}

# give it the dish surface accuracy in mm
self.getetaa := function(eta0=0.80, surfacemm=25) {
  wider private;
  lambdamm := 300.0/private.freq;
  sigmamm  := surfacemm/1000;
  private.etaa := eta0 * e^(- ( 4*pi*sigmamm / lambdamm)^2 );
  
  return private.etaa;
}

self.usejanskies := function() {
	wider private;
	private.sensunits := 'Jy';
# 	print 'Sensitivity Units are Jy';
	return T;
}
self.usekelvins := function() {
	wider private;
	private.sensunits := 'K ';
# 	print 'Sensitivity Units are K';
	return T;
}


# just set dt (integration time)
self.setdt := function(dt=1) {
	wider private;
	private.dt := dt;
	return T;
}
# set the details of the observation; 
#        freqs in GHz,  time in seconds, positions in degrees
self.setobservation := function(npol=2, bandwidth=8, freq=300, dt=1, rxspec=T) {
	wider private;
	private.npol := npol;
	private.bandwidth := bandwidth;
	private.freq := freq;
	private.dt := dt;
	private.rxspec := rxspec;
	return self.isinband();
}


self.isinband := function() {
	wider private;
	whichband := self.whichband( private.freq );
	if (whichband == 0) {
	   note(paste('ERROR: freq ', private.freq, ' is out of band'));
	   private.band := F;
	   private.trx  := F;
#	   fail;
	} else {
	   private.band := whichband;
	   if (private.rxspec) {
	     private.trx  := private.bandtrxspec[ whichband ];
	   } else {
	     private.trx  := private.bandtrxgoal[ whichband ];
	   }
	   return T;
	}
}

# returns the band number, or 0 if freq is not in any band
self.whichband := function(freq) {
	wider private;
# strategy:  pick lowest band this freq is in
	whichband := 0;
	for (i in [len(private.bands):1]) {
	   if (freq >= private.bandlow[i] && freq <= private.bandhigh[i]) {
	      whichband := i;
	   }
	}
	return whichband;
}


self.setinstrument := function(nants=1, diam=12, baseline=12) {
	wider private;
	private.nants := nants;
	private.diam := diam;
	private.baseline := baseline;
	return T;
}

# give either tau225 OR pwv (OR neither works if you use settaufreq())
self.setatmos := function(tatmos=250, tamb=270, tau225=F, pwv=F, airmass=1) {
	wider private;

	private.tatmos := tatmos;
	private.tamb := tamb;
	if (pwv!=F) {
	  private.pwv := pwv; 
	  private.tau225 := self.gettau(225, pwv)
	} else if (tau225!=F){
	  private.tau225 := tau225;
	  private.pwv := self.getpwv(tau225);
	}
	if (private.freq != F) {
	   private.taufreq := self.gettaufreq(tau225, private.freq);
 	}
	private.airmass := airmass;
	if (! private.silent) {
	   print 'Using pwv =', private.pwv, 'and tau225 =', private.tau225;
	}
	return T;
}

self.settaufreq := function(taufreq=0.1, freq=300) {
	wider private;

	private.freq := freq;
	private.taufreq := taufreq;
	wet := self.getwet(freq);
	dry := self.getdry(freq);
	private.pwv := (taufreq - dry)/wet;
	private.tau225 := self.gettau(freq=225, pwv=private.pwv);
	if (! private.silent) {
 	   print 'Using pwv =', private.pwv, 'and tau225 =', private.tau225;
	}
	return T;
}


self.plottransmission := function(pwv=1) {
	return (private.almatau.plottrans(pwv));
}


self.summary := function() {
  note('Summary of sensitivity parameters:');
  note('---------------------------------------');
  note(paste('  private.nants   = ',   private.nants ));
  note(paste('  private.diam    = ',   private.diam  ));
  note(paste('  private.baselin = ',   private.baseline ));
  note(paste('  private.freq    = ',   private.freq ));
  note(paste('  private.band    = ',   private.band ));
  note(paste('  private.trx     = ',   private.trx ));
  note('---------------------------------------');
  note(paste('  private.pwv     = ',   private.pwv ));
  note(paste('  private.tau225  = ',   private.tau225 ));
  note(paste('  private.taufreq = ',   private.taufreq ));
  note(paste('  private.tatmos  = ',   private.tatmos ));
  note(paste('  private.tamb    = ',   private.tamb   ));
  note(paste('  private.tcmb    = ',   private.tcmb   ));
  note(paste('  private.airmass = ',   private.airmass ));
  note(paste('  private.gaincor = ',   private.gaincorrect ));
  note('---------------------------------------');
  note(paste('  private.etaa    = ',   private.etaa ));
  note(paste('  private.etal    = ',   private.etal ));
  note(paste('  private.etaq    = ',   private.etaq ));
  note('---------------------------------------');
  note(paste('  private.npol    = ',   private.npol ));
  note(paste('  private.bandwdt = ',   private.bandwidth ));
  note(paste('  private.freq    = ',   private.freq ));
  note(paste('  private.dt      = ',   private.dt ));
  note('---------------------------------------');
  note(paste('  private.h       = ',   private.h ));
  note(paste('  private.k       = ',   private.k ));
  note(paste('  private.units   = ',   private.sensunits ));
  note('---------------------------------------');

}


self.calctsys := function() {

  private.verify();

  if (private.gaincorrect) {
    tterm1 := private.trx* e^( private.taufreq * private.airmass);
    tterm2 := private.etal * 
		private.tatmos * (e^(private.taufreq * private.airmass) - 1);
    tterm3 := (1 - private.etal)*private.tamb* 
		e^(private.taufreq * private.airmass);
    tterm4 := private.tcmb;
  } else {
    tterm1 := private.trx;
    tterm2 := private.etal * 
		private.tatmos * (1-e^(-private.taufreq * private.airmass));
    tterm3 := (1 - private.etal)*private.tamb;
    tterm4 := private.tcmb * e^(-private.taufreq * private.airmass);
  }

  tsys := tterm1 + tterm2 + tterm3 + tterm4;

  return tsys;
}




#  This is a high level sensitivity calculation that
#  just takes the atmospheric conditions and the band
#
self.calcsensitivity1 := function(mode='obs')  {
  
  # ds in Jy or K (switch using self.usekelvins() or self.usejanskies()
  ds := self.calcsensitivity0();
  
  # do the Nant thing
  if (mode == 'obs') {
    if (private.nants >= 2) {
      ds := ds / sqrt(private.nants*(private.nants-1)/2);
    } else {                 #  assume 1 antenna
      ds := ds/sqrt(2);   
      # Single Dish:  this does not include atmospheric switching
      # usually, SD has sqrt(2) up top: half time ON, 
      # and diff of two noisy numbers
    }
  } else if (mode == 'gain') { 
    if (private.nants >= 3) {
      ds := ds / sqrt(private.nants-2);
    } else {
      note(paste('Unknown mode', mode, 'in calcsensitivity1()'));
      return F;
    }
  }
  return (paste( ds, private.sensunits));
}



#  Calculate the point source noise for one baseline;
#  Low level routine with no assumptions about how things vary with
#  frequency
#
self.calcsensitivity0 := function() {

  private.verify();

  npol := private.npol;
  dvhz := private.bandwidth * 1e+9;

  tsys := self.calctsys();

  area := (private.diam)^2 * pi / 4;

  ds := ( 2^(0.5) * private.k * private.jy * tsys ) /
	( private.etaq * private.etaa * area *
	  sqrt(private.npol * dvhz * private.dt) );
  if (private.sensunits == 'K ') {
	ds := ds /(2 * private.k * private.jy * 1.316 / 
		(private.baseline)^2 ); 	    
  }
#  note: 1.316 is from 1.13309 * (theta/lambda)^2
#      and FIT theta == 1.16 lambda / D (for gaussian single
#      dish beam 
  return ds;
}

#  This is a high level sensitivity calculation;
#  It returns the integration time required to achieve the
#  desired sensitivity

self.calcinttime1 := function(sigma=0.01, mode='obs')  {

  # mydt in [s], for one baseline
  mydt := self.calcinttime0(sigma);
  
  # do the Nant thing
  if (mode == 'obs') {
    if (private.nants >= 2) {
      mydt := mydt / (private.nants*(private.nants-1)/2);
    } else {
      mydt := mydt;   #  this does not include atmospheric switching
    }
  } else if (mode == 'gain') { 
    if (private.nants >= 3) {
      mydt := mydt / (private.nants-2);
    } else {
      note(paste('Unknown mode', mode, 'in calcinttime1()'));
      return F;
    }
  }
  return mydt;
}

#  Calculate the time for one baseline to reach a given sigma
#  Low level routine with no assumptions about how things vary with
#  frequency
#
  self.calcinttime0 := function(sigma) {
	wider private;
	private.verify();

	savetime := private.dt;
	private.dt := 1.0;
	aa := self.calcsensitivity0();
	private.dt := savetime;
	mydt := (aa / sigma)^2
	return mydt;
  }



  self.gettau := function(freq=225, pwv=1) {
     return (private.almatau.gettau(freq, pwv));
  }
  self.getpwv := function(tau225=1) {
     return (private.almatau.getpwv(tau225));
  }
  self.getwet := function(freq=225) {
     return (private.almatau.getwet(freq));
  }
  self.getdry := function(freq=225) {
     return (private.almatau.getdry(freq));
  }
  self.gettaufreq := function(tau225, freq) {
	pwv := self.getpwv(tau225);
	return (self.gettau(freq, pwv));
  }

##################################################################################
######### Private Functions: reservations must be made 2 weeks in advance ########
##################################################################################

  private.verify := function() {
    wider private;
    if (private.freq==F) {
      print 'Common, meet me half way: what FREQUENCY do you want?';
      fail;
    }
    private.taufreq := self.gettaufreq(private.tau225, private.freq);

    return T;
  }


# end of constructor
}



testalmasensitivity := function(units='Jy') {


  include 'almasensitivity.g'
  mysens := almasensitivity();
  if (units=='Jy') {
     mysens.usejanskies();
  } else {
     mysens.usekelvins();
  }

  mysens.setinstrument(nants=64, diam=12, baseline=12);
  mysens.setatmos(tau225=0.05);
  mysens.setatmos(pwv=1.3446, airmass=1.3);

  freq := 35;  tau := 0.017;
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60);
  print 'band is: ', mysens.whichband(freq)
  mysens.settaufreq(freq=freq, taufreq=tau); 
  tsys := mysens.calctsys();
  sens := mysens.calcsensitivity1('obs');
  print 'FREQ = ', freq, '   tau = ', tau, '  Tsys = ', tsys, '  sigma = ', sens;

  freq := 110;  tau := 0.049;
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60);
  print 'band is: ', mysens.whichband(freq)
  mysens.settaufreq(freq=freq, taufreq=tau); 
  tsys := mysens.calctsys();
  sens := mysens.calcsensitivity1('obs');
  print 'FREQ = ', freq, '   tau = ', tau, '  Tsys = ', tsys, '  sigma = ', sens;

  freq := 230;  tau := 0.078;
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60);
  print 'band is: ', mysens.whichband(freq)
  mysens.settaufreq(freq=freq, taufreq=tau); 
  tsys := mysens.calctsys();
  sens := mysens.calcsensitivity1('obs');
  print 'FREQ = ', freq, '   tau = ', tau, '  Tsys = ', tsys, '  sigma = ', sens;

  freq := 345;  tau := 0.276;
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60);
  print 'band is: ', mysens.whichband(freq)
  mysens.settaufreq(freq=freq, taufreq=tau); 
  tsys := mysens.calctsys();
  sens := mysens.calcsensitivity1('obs');
  print 'FREQ = ', freq, '   tau = ', tau, '  Tsys = ', tsys, '  sigma = ', sens;

  freq := 409;  tau := 0.544;
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60);
  print 'band is: ', mysens.whichband(freq)
  mysens.settaufreq(freq=freq, taufreq=tau); 
  tsys := mysens.calctsys();
  sens := mysens.calcsensitivity1('obs');
  print 'FREQ = ', freq, '   tau = ', tau, '  Tsys = ', tsys, '  sigma = ', sens;

  freq := 675;  tau := 1.789;
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60);
  print 'band is: ', mysens.whichband(freq)
  mysens.settaufreq(freq=freq, taufreq=tau); 
  tsys := mysens.calctsys();
  sens := mysens.calcsensitivity1('obs');
  print 'FREQ = ', freq, '   tau = ', tau, '  Tsys = ', tsys, '  sigma = ', sens;

  freq := 860;  tau := 1.601; 
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60, rxspec=T);
  print 'band is: ', mysens.whichband(freq)
  mysens.settaufreq(freq=freq, taufreq=tau); 
  tsys := mysens.calctsys();
  sens := mysens.calcsensitivity1('obs');
  print 'FREQ = ', freq, '   tau = ', tau, '  Tsys = ', tsys, '  sigma = ', sens;

  mysens.summary();

  mysens.done();

}