#
#
#   gets data to make a plot of residuals vs switching rate
#
#
include 'sdatmosim.g'
include 'almatau.g'

iband := 6
el    := 60
gflucs := 0.0

tau225good :=
[0.328999996,0.175999999,0.116999999,0.0839999989,0.0390000008,0.057,0.0450000018,0.0359999985,0.0309999995,0.023]  
phagood := 
[5.19000006,4.11999989,3.26999998,2.11999989,5.09000015,1.97000003,1.71500003,1.44000006,1.00999999,0.689999998]
tau225bad  :=
[0.138999999,0.129999995,0.111000001,0.0890000015,0.0549999997,0.057,0.0480000004,0.0399999991,0.0340000018,0.0270000007]  
phabad  :=
[8.29500008,5.73000002,4.2249999,3.1500001,4.31500006,2.6099999,2.05500007,1.58000004,1.23500001,0.879999995]
freq := [43, 80, 90, 145, 190, 230, 345, 500, 680, 880];

  mytautool := almatau();
  pwv := mytautool.getpwv(tau225good[iband]);
  mysim := sdatmosim()
  mysim.readaips('ATMO.DOUBLE')
  mysim.setup(azimuth=225, elevation=el, pwv=pwv, stiphase=0, 
		freq=freq[iband], gainfluctuations=0.0)


dtvec := [0.01, 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4];
residnoise := 0 * dtvec;
residatmo  := 0 * dtvec;
residgain3 := 0 * dtvec;
residgain4 := 0 * dtvec;


for (ii in ind(dtvec)) {
  nb := as_integer( 10 * .3 / dtvec[ii] )
  print ii, nb
  mysim.setupbs2(mythrow=1, dt=dtvec[ii], tset=0.005, nsmallcycles=10, nbigcycles=nb)
  results := mysim.bs2(F, T);
  residnoise[ii] := results.rms
}  


mysim.setup(azimuth=225, elevation=el, pwv=pwv, stiphase=phagood[iband], 
		freq=freq[iband], gainfluctuations=0.0)

for (ii in ind(dtvec)) {
  nb := as_integer( 10 * .3 / dtvec[ii] )
  print ii, nb
  mysim.setupbs2(mythrow=1, dt=dtvec[ii], tset=0.005, nsmallcycles=10, nbigcycles=nb)
  results := mysim.bs2(T, T);
  residatmo[ii] := results.rms
}  



mysim.setup(azimuth=225, elevation=el, pwv=pwv, stiphase=phagood[iband], 
		freq=freq[iband], gainfluctuations=0.001 )

for (ii in ind(dtvec)) {
  nb := as_integer( 10 * .3 / dtvec[ii] )
  print ii, nb
  mysim.setupbs2(mythrow=1, dt=dtvec[ii], tset=0.005, nsmallcycles=10, nbigcycles=nb)
  results := mysim.bs2(F, T);
  residgain3[ii] := results.rms
}  


mysim.setup(azimuth=225, elevation=el, pwv=pwv, stiphase=phagood[iband], 
		freq=freq[iband], gainfluctuations=0.0001 )

for (ii in ind(dtvec)) {
  nb := as_integer( 10 * .3 / dtvec[ii] )
  print ii, nb
  mysim.setupbs2(mythrow=1, dt=dtvec[ii], tset=0.005, nsmallcycles=10, nbigcycles=nb)
  results := mysim.bs2(F, T);
  residgain4[ii] := results.rms
}  


for (ii in ind(dtvec)) {
  print  ii,  residnoise[ii],  residatmo[ii],  residgain3[ii],  residgain4[ii]
}

mytautool.done();
mysim.done();