#
#
#  Strategy:  drive the fastswitchsim2
#      1) get the atmospheric parameters correct
#      2) drive ~1000 monte carlo sims of the atmosphere
#          for each frequency & elevation
#      3) capture results and take medians
#      4) present results (yeah, RIGHT!)

include 'sourcecountsim.g'
include 'fastswitchsim2.g'





#  flat 8% per band;
#  tau == tau225
#  pha == rms phase @ 300m, 11.2 GHz
#
#fr band -median--  -- 3/4 ---
#       tau   pha   tau   pha
#
# 43  1 0.22  6.71  0.116 8.28
# 80  2 0.129 5.44  0.108 6.125
# 90  3 0.129 3.635 0.119 4.315
#145  4 0.089 2.955 0.087 3.58
#190  5 0.06  3.21  0.065 3.295
#230  6 0.056 2.175 0.057 2.655
#345  7 0.045 1.84  0.047 2.015
#500  8 0.037 1.615 0.04  1.625
#680  9 0.032 1.155 0.034 1.255
#880 10 0.025 0.79  0.027 0.89
#
#unused 0.101 11.44 0.092 12.775

#  TEST
#
# driveone(targetfreq= 230, calfreq= 230, el=60, tau225=0.045, rms300=1.84, 
#    outfile='AAA.TEST3', nsim=5)


driveallc := function() {

  driveone(targetfreq=850, calfreq= 90, el=60, tau225=0.032, kinkalpha=0.5, rms300=0.790, outfile='AAC.850.60')
  driveone(targetfreq=650, calfreq= 90, el=60, tau225=0.037, kinkalpha=0.5, rms300=1.155, outfile='AAC.650.60')
  driveone(targetfreq=490, calfreq= 90, el=60, tau225=0.045, kinkalpha=0.5, rms300=1.615, outfile='AAC.490.60')
  driveone(targetfreq=345, calfreq= 90, el=60, tau225=0.045, kinkalpha=0.5, rms300=1.840, outfile='AAC.345.60')
  driveone(targetfreq=230, calfreq= 90, el=60, tau225=0.056, kinkalpha=0.5, rms300=2.175, outfile='AAC.225.60')
  driveone(targetfreq=230, calfreq=230, el=60, tau225=0.056, kinkalpha=0.5, rms300=2.175, outfile='AAC.225B.60')
  driveone(targetfreq=190, calfreq= 90, el=60, tau225=0.060, kinkalpha=0.5, rms300=3.210, outfile='AAC.190.60')
  driveone(targetfreq=190, calfreq=190, el=60, tau225=0.060, kinkalpha=0.5, rms300=3.210, outfile='AAC.190B.60')
  driveone(targetfreq=145, calfreq= 90, el=60, tau225=0.089, kinkalpha=0.5, rms300=2.955, outfile='AAC.145.60')
  driveone(targetfreq=145, calfreq=145, el=60, tau225=0.089, kinkalpha=0.5, rms300=2.955, outfile='AAC.145B.60')
  driveone(targetfreq= 90, calfreq= 90, el=60, tau225=0.129, kinkalpha=0.5, rms300=3.635, outfile='AAC.090.60')
  driveone(targetfreq= 70, calfreq= 70, el=60, tau225=0.129, kinkalpha=0.5, rms300=5.440, outfile='AAC.070.60')
  driveone(targetfreq= 35, calfreq= 35, el=60, tau225=0.220, kinkalpha=0.5, rms300=6.710, outfile='AAC.035.60')
  driveone(targetfreq=850, calfreq=850, el=60, tau225=0.032, kinkalpha=0.5, rms300=0.790, outfile='AAC.850B.60')
  driveone(targetfreq=650, calfreq=650, el=60, tau225=0.037, kinkalpha=0.5, rms300=1.155, outfile='AAC.650B.60')
  driveone(targetfreq=490, calfreq=490, el=60, tau225=0.045, kinkalpha=0.5, rms300=1.615, outfile='AAC.490B.60')
  driveone(targetfreq=345, calfreq=345, el=60, tau225=0.045, kinkalpha=0.5, rms300=1.840, outfile='AAC.345B.60')

  return T;

}

leftoversc := function() {

  driveone(targetfreq=145, calfreq= 90, el=60, tau225=0.089, kinkalpha=0.5, rms300=2.955, outfile='AAC.145.60')
  driveone(targetfreq=145, calfreq=145, el=60, tau225=0.089, kinkalpha=0.5, rms300=2.955, outfile='AAC.145B.60')
  driveone(targetfreq= 90, calfreq= 90, el=60, tau225=0.129, kinkalpha=0.5, rms300=3.635, outfile='AAC.090.60')
  driveone(targetfreq= 70, calfreq= 70, el=60, tau225=0.129, kinkalpha=0.5, rms300=5.440, outfile='AAC.070.60')
  driveone(targetfreq= 35, calfreq= 35, el=60, tau225=0.220, kinkalpha=0.5, rms300=6.710, outfile='AAC.035.60')
  driveone(targetfreq=850, calfreq=850, el=60, tau225=0.032, kinkalpha=0.5, rms300=0.790, outfile='AAC.850B.60')
  driveone(targetfreq=650, calfreq=650, el=60, tau225=0.037, kinkalpha=0.5, rms300=1.155, outfile='AAC.650B.60')
  driveone(targetfreq=490, calfreq=490, el=60, tau225=0.045, kinkalpha=0.5, rms300=1.615, outfile='AAC.490B.60')
  driveone(targetfreq=345, calfreq=345, el=60, tau225=0.045, kinkalpha=0.5, rms300=1.840, outfile='AAC.345B.60')

  return T;

}


drivealld := function() {

  driveone(targetfreq=850, calfreq= 90, el=60, tau225=0.032, kinkalpha=0.0, rms300=0.790, outfile='AAD.850.60')
  driveone(targetfreq=650, calfreq= 90, el=60, tau225=0.037, kinkalpha=0.0, rms300=1.155, outfile='AAD.650.60')
  driveone(targetfreq=490, calfreq= 90, el=60, tau225=0.045, kinkalpha=0.0, rms300=1.615, outfile='AAD.490.60')
  driveone(targetfreq=345, calfreq= 90, el=60, tau225=0.045, kinkalpha=0.0, rms300=1.840, outfile='AAD.345.60')
  driveone(targetfreq=230, calfreq= 90, el=60, tau225=0.056, kinkalpha=0.0, rms300=2.175, outfile='AAD.225.60')
  driveone(targetfreq=230, calfreq=230, el=60, tau225=0.056, kinkalpha=0.0, rms300=2.175, outfile='AAD.225B.60')
  driveone(targetfreq=190, calfreq= 90, el=60, tau225=0.060, kinkalpha=0.0, rms300=3.210, outfile='AAD.190.60')
  driveone(targetfreq=190, calfreq=190, el=60, tau225=0.060, kinkalpha=0.0, rms300=3.210, outfile='AAD.190B.60')
  driveone(targetfreq=145, calfreq= 90, el=60, tau225=0.089, kinkalpha=0.0, rms300=2.955, outfile='AAD.145.60')
  driveone(targetfreq=145, calfreq=145, el=60, tau225=0.089, kinkalpha=0.0, rms300=2.955, outfile='AAD.145B.60')
  driveone(targetfreq= 90, calfreq= 90, el=60, tau225=0.129, kinkalpha=0.0, rms300=3.635, outfile='AAD.090.60')
  driveone(targetfreq= 70, calfreq= 70, el=60, tau225=0.129, kinkalpha=0.0, rms300=5.440, outfile='AAD.070.60')
  driveone(targetfreq= 35, calfreq= 35, el=60, tau225=0.220, kinkalpha=0.0, rms300=6.710, outfile='AAD.035.60')
  driveone(targetfreq=850, calfreq=850, el=60, tau225=0.032, kinkalpha=0.0, rms300=0.790, outfile='AAD.850B.60')
  driveone(targetfreq=650, calfreq=650, el=60, tau225=0.037, kinkalpha=0.0, rms300=1.155, outfile='AAD.650B.60')
  driveone(targetfreq=490, calfreq=490, el=60, tau225=0.045, kinkalpha=0.0, rms300=1.615, outfile='AAD.490B.60')
  driveone(targetfreq=345, calfreq=345, el=60, tau225=0.045, kinkalpha=0.0, rms300=1.840, outfile='AAD.345B.60')

  return T;

}



makefields := function(nsim=1000) {

  for (i in [1:nsim]) {
    mycounts := sourcecountsim();
    mycounts.setsourcecounts(smin=0.005, fieldsize=40);  
    mycounts.setupcounts();
    mycounts.readalpha();
    sourcelist :=  mycounts.simulatesources();
    mycounts.done();

    filename := getfilename(i);
    record2table( sourcelist, filename, 'sourcelist');
    sourcelist := F;
  }
  return T;

}


one_for_larry := function(outfile='LARRY.SUM', distfile='LARRY.DIST') {

   driveone(targetfreq=937, calfreq=90, el=60, tau225=0.0296
                     kinkalpha=0.5, rms300=0.927, alpha=0.58, nsim=1000,
		     outfile=spaste(outfile,'.937'), 
		     distfile=spaste(distfile, '.937'));

   driveone(targetfreq=875, calfreq=90, el=60, tau225=0.0296
                     kinkalpha=0.5, rms300=0.927, alpha=0.58, nsim=1000,
		     outfile=spaste(outfile,'.875'), 
		     distfile=spaste(distfile, '.875'));

   print 'Look in these files:'
   print spaste(outfile,'.937'), spaste(distfile, '.937'),
         spaste(outfile,'.875'), spaste(distfile, '.875');
   return T;
}

five_for_sensitivity := function(freq=875) {

  mysens := almasensitivity();
  mysens.usejanskies();

  mysens.setinstrument(nants=64, diam=12, baseline=12);
  mysens.setobservation(npol=2, bandwidth=8, freq=freq, dt=60);
  print " tau   tsys   sens ";
  for (tau in [0.02875,0.0296,0.0304,0.0326,0.0348]) {
    mysens.setatmos(tau225=tau, airmass=1/cos(30/57.29));
    tsys := mysens.calctsys();
    sens := mysens.calcsensitivity1('obs');
    print tau, tsys, sens;
  }
  mysens.done();
  return T
}


five_for_larry := function(outfile='LARRY.SUM', distfile='LARRY.DIST') {

   driveone(targetfreq=937, calfreq=90, el=60, tau225=0.02875,
                     kinkalpha=0.5, rms300=0.985, alpha=0.58, nsim=1000,
		     outfile=spaste(outfile,'.1.937'), 
		     distfile=spaste(distfile, '.1.937'));

   driveone(targetfreq=937, calfreq=90, el=60, tau225=0.0296,
                     kinkalpha=0.5, rms300=0.927, alpha=0.58, nsim=1000,
		     outfile=spaste(outfile,'.2.937'), 
		     distfile=spaste(distfile, '.2.937'));

   driveone(targetfreq=937, calfreq=90, el=60, tau225=0.0304,
                     kinkalpha=0.5, rms300=0.869, alpha=0.58, nsim=1000,
		     outfile=spaste(outfile,'.3.937'), 
		     distfile=spaste(distfile, '.3.937'));

   driveone(targetfreq=937, calfreq=90, el=60, tau225=0.0326,
                     kinkalpha=0.5, rms300=0.811, alpha=0.58, nsim=1000,
		     outfile=spaste(outfile,'.4.937'), 
		     distfile=spaste(distfile, '.4.937'));

   driveone(targetfreq=937, calfreq=90, el=60, tau225=0.0348,
                     kinkalpha=0.5, rms300=0.753, alpha=0.58, nsim=1000,
		     outfile=spaste(outfile,'.5.937'), 
		     distfile=spaste(distfile, '.5.937'));


   print 'Look in these files:'
   print spaste(outfile,'.1.937'), spaste(outfile,'.2.937'), 
	spaste(outfile,'.3.937'), spaste(outfile,'.4.937'), 
	spaste(outfile,'.5.937');
   print spaste(distfile,'.1.937'), spaste(distfile,'.2.937'), 
	spaste(distfile,'.3.937'), spaste(distfile,'.4.937'), 
	spaste(distfile,'.5.937'); 

   return T;
}

#
driveone := function(targetfreq=650, calfreq=90, el=60, tau225=0.05,
                     kinkalpha=0.5, rms300=1.0, alpha=0.58, nsim=1000,
		     outfile='TEST.OUT', distfile=' ') {

# what parameters to store
  
  all_phase    := [1:nsim]
  all_eff0     := [1:nsim]
  all_eff1     := [1:nsim]
  all_eff2     := [1:nsim]
  cb_eff0      := [1:nsim]
  cb_eff1      := [1:nsim]
  cb_eff2      := [1:nsim]
  fs_eff0      := [1:nsim]
  fs_eff1      := [1:nsim]
  fs_eff2      := [1:nsim]
  fs_calflux   := [1:nsim]
  fs_caldist   := [1:nsim]
  fs_tslew     := [1:nsim]
  fs_tcal      := [1:nsim]
  fs_tcycle    := [1:nsim]
  fs_ttarget   := [1:nsim]
  fs_toverhead := [1:nsim]
  fs_phase     := [1:nsim]
  fs_atmophase := [1:nsim]
  fs_gainerror := [1:nsim]  
  fs_cb_samecal:= 0 * [1:nsim]  
  cb_calflux   := [1:nsim]
  cb_tarflux   := [1:nsim]
  cb_caldist   := [1:nsim]
  cb_tslew     := [1:nsim]
  cb_tcal      := [1:nsim]
  cb_tcycle    := [1:nsim]
  cb_ttarget   := [1:nsim]
  cb_toverhead := [1:nsim]
  cb_phase     := [1:nsim]
  cb_atmophase := [1:nsim]
  cb_gainerror := [1:nsim]  

  for (i in [1:nsim]) {

    filename := getfilename(i);
    if ( ! dos.fileexists(filename) )  {
       print 'HEY, source list file', filename, 'not found';
       print 'Better run makefields(nsim=', nsim, ')';
       return F;
    }

    sourcelist :=  table2record( filename, type='sourcelist' )
    if (! is_record( sourcelist) ) {
       print 'HEY, sourcelist did not read correctly from', filename;
       return F;
    }
    
    sim := fastswitchsim2()
    sim.setatmos(rms=rms300, alpha=alpha, tau225=tau225)
    sim.setkinkalpha( kinkalpha );
        
    sdat := sim.optimizeswitching(sourcelist=sourcelist,
    calfreq=calfreq, targetfreq=targetfreq, tchangefreq=1.0, 
    el=el, velocity=12, crosscycletime=300) 

    sim.done();
    sourcelist := F;

    all_phase[i]    := sdat.all.residphase;
    all_eff0[i]     := sdat.all.eff0;
    all_eff1[i]     := sdat.all.eff1;
    all_eff2[i]     := sdat.all.eff2;
    cb_eff0[i]      := sdat.cb.eff0;
    cb_eff1[i]      := sdat.cb.eff1;
    cb_eff2[i]      := sdat.cb.eff2;
    fs_eff0[i]      := sdat.fs.eff0;
    fs_eff1[i]      := sdat.fs.eff1;
    fs_eff2[i]      := sdat.fs.eff2;
    fs_calflux[i]   := sdat.fs.calflux;
    fs_caldist[i]   := sdat.fs.caldistance;
    fs_tslew[i]     := sdat.fs.tslew;
    fs_tcal[i]      := sdat.fs.tcal;
    fs_tcycle[i]    := sdat.fs.tcycle;
    fs_ttarget[i]   := sdat.fs.ttarget;
    fs_toverhead[i] := sdat.fs.toverhead;
    fs_phase[i]     := sdat.fs.residphase;
    fs_atmophase[i] := sdat.fs.atmophase;
    fs_gainerror[i] := sdat.fs.gainerror;

    cb_calflux[i]   := sdat.cb.calfluxcal;
    cb_tarflux[i]   := sdat.cb.calfluxtarget;
    cb_caldist[i]   := sdat.cb.caldistance;
    cb_tslew[i]     := sdat.cb.tslew;
    cb_tcal[i]      := sdat.cb.tcal;
    cb_tcycle[i]    := sdat.cb.tcycle;
    cb_ttarget[i]   := sdat.cb.ttarget;
    cb_toverhead[i] := sdat.cb.toverhead;
    cb_phase[i]     := sdat.cb.residphase;
    cb_atmophase[i] := sdat.cb.atmophase;
    cb_gainerror[i] := sdat.cb.gainerror;

    if ( sdat.fs.calflux == sdat.cb.calfluxcal ) {
       fs_cb_samecal[i]:= 1;
    }
  }
# end of Monte Carlo Sim

  out := open(paste( ">", outfile));

  write(out, 'targetfreq   ',      targetfreq, sep=' ');
  write(out);
  write(out, 'calfreq      ',      calfreq, sep=' ');
  write(out);
  write(out, 'el           ',      el, sep=' ');
  write(out);
  write(out, 'tau225       ',      tau225, sep=' ');
  write(out);
  write(out, 'rms300       ',      rms300, sep=' ');
  write(out);
  write(out, '###############################################');
  write(out);
  write(out, 'all_phase    ',      median(all_phase) , sep=' ');      
  write(out);
  write(out, 'all_eff0     ',      median(all_eff0) , sep=' ');       
  write(out);
  write(out, 'all_eff1     ',      median(all_eff1) , sep=' ');       
  write(out);
  write(out, 'all_eff2     ',      median(all_eff2) , sep=' ');       
  write(out);
  write(out, 'cb_eff0      ',      median(cb_eff0) , sep=' ');        
  write(out);
  write(out, 'cb_eff1      ',      median(cb_eff1) , sep=' ');        
  write(out);
  write(out, 'cb_eff2      ',      median(cb_eff2) , sep=' ');        
  write(out);
  write(out, 'fs_eff0      ',      median(fs_eff0) , sep=' ');        
  write(out);
  write(out, 'fs_eff1      ',      median(fs_eff1) , sep=' ');        
  write(out);
  write(out, 'fs_eff2      ',      median(fs_eff2) , sep=' ');        
  write(out);
  write(out, 'fs_calflux   ',      median(fs_calflux) , sep=' ');     
  write(out);
  write(out, 'fs_caldist   ',      median(fs_caldist) , sep=' ');     
  write(out);
  write(out, 'fs_tslew     ',      median(fs_tslew) , sep=' ');       
  write(out);
  write(out, 'fs_tcal      ',      median(fs_tcal) , sep=' ');        
  write(out);
  write(out, 'fs_tcycle    ',      median(fs_tcycle) , sep=' ');      
  write(out);
  write(out, 'fs_ttarget   ',      median(fs_ttarget) , sep=' ');     
  write(out);
  write(out, 'fs_toverhead ',      median(fs_toverhead) , sep=' ');   
  write(out);
  write(out, 'fs_phase     ',      median(fs_phase) , sep=' ');       
  write(out);
  write(out, 'fs_atmophase ',      median(fs_atmophase) , sep=' ');       
  write(out);
  write(out, 'fs_gainerror ',      median(fs_gainerror) , sep=' ');      
  write(out);
  write(out, 'cb_calflux   ',      median(cb_calflux) , sep=' ');     
  write(out);
  write(out, 'cb_tarflux   ',      median(cb_tarflux) , sep=' ');     
  write(out);
  write(out, 'cb_caldist   ',      median(cb_caldist) , sep=' ');     
  write(out);
  write(out, 'cb_tslew     ',      median(cb_tslew) , sep=' ');       
  write(out);
  write(out, 'cb_tcal      ',      median(cb_tcal) , sep=' ');        
  write(out);
  write(out, 'cb_tcycle    ',      median(cb_tcycle) , sep=' ');      
  write(out);
  write(out, 'cb_ttarget   ',      median(cb_ttarget) , sep=' ');     
  write(out);
  write(out, 'cb_toverhead ',      median(cb_toverhead) , sep=' ');   
  write(out);
  write(out, 'cb_phase     ',      median(cb_phase) , sep=' ');       
  write(out);
  write(out, 'cb_atmophase ',      median(cb_atmophase) , sep=' ');       
  write(out);
  write(out, 'cb_gainerror ',      median(cb_gainerror) , sep=' ');      
  write(out);
  write(out, 'Same cal     ', sum(fs_cb_samecal)/len(fs_cb_samecal), sep=' ');
  write(out);
  out := F

  if (distfile != ' ') {
#  write uot distributions

       out := open(paste( ">", distfile));
       write(out, '# phase  eff0  eff1  eff2  --FAST-SWITCHING--    --CROSS-BAND-CAL--');
       write(out);
       write(out, '#                          phase calflx caldist  phase calflx tarflx caldist');
       write(out);
       write(out, '#  1       2      3     4     5     6     7      8      9      10      11');
       write(out);
       for (i in ind(all_phase)) {
          write(out, huh(all_phase[i],1), huh(all_eff0[i],3), huh(all_eff1[i],3), 
                     huh(all_eff2[i],3), huh(fs_phase[i],1), huh(fs_calflux[i],3), 
                     huh(fs_caldist[i],2), huh(cb_phase[i],1), huh(cb_calflux[i],3), 
                     huh(cb_tarflux[i],3), huh(cb_caldist[i],2), sep=' ');
          write(out);
       }
       out := F;
  }

  return T;
}

huh := function(anum, nn=3) {
  aa := 10^nn
  return( as_integer( aa * anum ) / aa );
}

getfilename := function( i ) {

   dirname := 'TEMP_SOURCELISTS';
   if ( ! dos.fileexists(dirname) ) {
      shell (paste( ' mkdir ', dirname ));
   }

   filename := spaste( dirname, '/', 'field.', i)
   return filename;
}

# write record to a table
record2table := function( ref myrecord=[=], 
                          tablename='', 
                          type='unknown' ){

    myrecord.type := type;
    scd:=tablecreatescalarcoldesc("dummy",[=]);
    td:=tablecreatedesc(scd);
    tab:=table(tablename, tabledesc=td, nrow=0)
    tab.putkeyword("datarecord", myrecord);
    tab.done();
    return T;	
}
  
# read in a record from a table
#
table2record := function( tablename='', type='') {

    tab:=table(tablename);
    myrecord := tab.getkeyword("datarecord");
    if (is_fail(myrecord)) {
      note(paste("Failed to read table", tablename));
      fail;
    }
    if (type != '') {
      if (myrecord.type != type) {
	note(paste("table", tablename, "is not of correct type", type));
	fail;
      }
    }	
    tab.done();
    
    return ref myrecord;
}