# fastswitchsim2.g: definition of fast switching simulator # ## 2003dec13 LAMA Memo 803 writeup version # # # sim := fastswitchsim2() # sim.setatmos(rms=1.0, alpha=0.6, tau225=0.05) # setkinkalpha( alpha=0.5 ) # # # pragma include once; include 'quanta.g' include 'measures.g' include 'pgplotter.g' include 'stidata.g' include 'almasensitivity.g' include 'sourcecountsim.g' #const fastswitchsim2 := subsequence () { private := [=]; private.mysens := almasensitivity(); private.mysens.usejanskies(); private.mysens.setinstrument(nants=64, diam=12, baseline=12); private.mysens.silent(); private.stidata := stidata(); private.gamma0 := 0.93; # fraction of structure function for rms private.gamma1 := 0.50; # fraction os structure function for ave phase error private.h := 500; # height of turbulent layer [m] private.tau225 := 0.05; # opacity at 225 GHz private.mypg := F; private.kinkalpha := 0.5; ############################################################# # # Public functions # ############################################################# # # for debugging purposes # self.getprivate := function() { return ref private; } self.done := function() { wider self, private; if (is_tool(private.mypg) ) { private.mypg.done(); } private.mysens.done(); val self := F; val public := F; return T; } # # change ths spectral index for sources above 90 GHz # self.setkinkalpha := function( alpha=0.5 ) { wider private; private.kinkalpha := alpha; return T; } self.setatmos := function(rms=1.0, alpha=0.6, tau225=0.05) { wider private; private.tau225 := tau225; return private.stidata.setdata(rms=rms, alpha=alpha); } # # given observing and atmospheric parameters, determine the # optimum cal strategy # # # self.optimizeswitching := function(sourcelist=[=], calfreq=90, targetfreq=345, tchangefreq=1.0, el=60, velocity=12, crosscycletime=300) { gainerror := 5.0; # degrees; initial guess; dvtrecord := self.dvt(sourcelist0=sourcelist, calfreq=calfreq, targetfreq=targetfreq, tchangefreq=tchangefreq, el=el, v=velocity, gainerror=gainerror, beta1=8.0); crossbandrecord := self.crossband(sourcelist=sourcelist, calfreq=calfreq, targetfreq=targetfreq, tchangefreq=tchangefreq, el=el, v=velocity, gainerror=gainerror, crosscycletime=crosscycletime); evalrecord := self.evaluate( dvtrecord, crossbandrecord, targetfreq, el, gainerror ); dvtrecord := self.dvt(sourcelist0=sourcelist, calfreq=calfreq, targetfreq=targetfreq, tchangefreq=tchangefreq, el=el, v=velocity, gainerror=evalrecord.fs.gainerror, beta1=evalrecord.beta); crossbandrecord := self.crossband(sourcelist=sourcelist, calfreq=calfreq, targetfreq=targetfreq, tchangefreq=tchangefreq, el=el, v=velocity, gainerror=crossbandrecord.cbgainerror, crosscycletime=crosscycletime); evalrecord := self.evaluate( dvtrecord, crossbandrecord, targetfreq, el, gainerror ); results := [=]; results.fs := [=]; # fast switching results.cb := [=]; # cross band results.all := [=]; # overall results.calfreq := calfreq; results.targetfreq := targetfreq; # fast switching results.fs.calflux := dvtrecord.s; results.fs.caldistance := sqrt( (dvtrecord.dx)^2 + (dvtrecord.dy)^2 ); results.fs.dvt := evalrecord.fs.optdvt; results.fs.tslew := dvtrecord.tslew; results.fs.tcal := dvtrecord.tcal; results.fs.toverhead := dvtrecord.toverhead; results.fs.tcycle := evalrecord.beta * dvtrecord.toverhead; results.fs.ttarget := results.fs.tcycle - dvtrecord.toverhead; results.fs.duty := (results.fs.tcycle - results.fs.toverhead) / results.fs.toverhead; results.fs.residphase := evalrecord.fs.optresid; results.fs.atmophase := evalrecord.fs.atmophase; # atmospheric contribution to phase results.fs.gainerror := evalrecord.fs.gainerror; # permitted gain error results.fs.eff1 := evalrecord.fs.eff1; # fs duty eff results.fs.eff2 := evalrecord.fs.eff2; # fs decor eff results.fs.eff0 := evalrecord.fs.eff0; # overall fs eff # cross band results.cb.calfluxtarget := crossbandrecord.calfluxtarget; # cb cal S @ targetfreq results.cb.calfluxcal := crossbandrecord.calfluxcal; # cb cal S @ calfreq results.cb.caldistance := crossbandrecord.caldistance results.cb.tslew := crossbandrecord.tslew; results.cb.nints := crossbandrecord.nints; results.cb.tcal := crossbandrecord.tcal; results.cb.toverhead := crossbandrecord.toverhead; results.cb.tcycle := crosscycletime; results.cb.ttarget := results.cb.tcycle - crossbandrecord.toverhead; results.cb.duty := crossbandrecord.duty; results.cb.residphase := crossbandrecord.optresid; results.cb.atmophase := crossbandrecord.atmophase; results.cb.atmophase := crossbandrecord.atmophase; results.cb.noisephase := crossbandrecord.noisephase; results.cb.gainerror := crossbandrecord.cbgainerror; results.cb.eff1 := crossbandrecord.eff1; # cb duty eff results.cb.eff2 := crossbandrecord.eff2; # cb decorr eff results.cb.eff0 := crossbandrecord.eff0; # overall cb eff # overall results.all.residphase := evalrecord.all.optresid; # overall rms phase results.all.eff1 := evalrecord.all.eff1; # overall duty eff results.all.eff2 := evalrecord.all.eff2; # overall decorr eff results.all.eff0 := evalrecord.all.eff0; # overall eff return ref results; } # # Figure out which sources will make a good cross band calibrator # and return info # self.crossband := function(sourcelist=[=], calfreq=90, targetfreq=345, tchangefreq=1, el=60, v=12, gainerror=5, crosscycletime=300) { results := [=]; if (targetfreq == calfreq) { # then we need NO crossband calibration results.needcrossband := F; results.eff0 := 1; results.eff1 := 1; results.duty := 1; results.eff2 := 1; results.optresid := 0; results.cbgainerror := 0; results.atmophase := 0; results.noisephase := 0; results.calfluxtarget := -1; results.calfluxcal := -1; results.caldistance := 0; results.nints := 0; results.toverhead := 0; results.tslew := 0; results.tcal := 0; } else { # first, select out just sources brighter than *fluxmin* mJy at the target freq results.needcrossband := T; fluxmin := 0.050; if (targetfreq > 90) { targetflux0 := sourcelist.s * (targetfreq/sourcelist.freq)^ (-sourcelist.alpha - private.kinkalpha); } else { targetflux0 := sourcelist.s * (targetfreq/sourcelist.freq)^(-sourcelist.alpha); } tforder0 := order(-targetflux0); tforder := tforder0[( targetflux0[tforder0] > fluxmin)]; if (len(tforder) == 0) { tforder[1] := tforder0[1]; } calflux := sourcelist.s[tforder]; targetflux := targetflux0[tforder]; x := sourcelist.x[tforder]; y := sourcelist.y[tforder]; nints0 := 0 * calflux; tgerr0 := 0 * calflux; toverhead0 := 0 * calflux; # our target is to have a residual gain error due to thermal noise of "gainerror" # Q: how long do we need to observe to get that? # sig1cal is the 1 second noise at the calibrator freq for (i in ind(tforder)) { private.mysens.setatmos(tau225=private.tau225, airmass=1/sin(el/(180/pi)) ); private.mysens.setobservation(npol=2, bandwidth=8, freq=calfreq, dt=1); sig1cal := 1.4142 * as_float( split( private.mysens.calcsensitivity1('gain'))[1] ); snr1cal := calflux[i] / sig1cal # remember: SNR of 10 in complex vector == 4.06 deg phase, rms # Question: how many cycles of 1s/1s/ do we need to get to "gainerror"? phase1cal := 4.06 * 10 / snr1cal; private.mysens.setobservation(npol=2, bandwidth=8, freq=targetfreq, dt=1); sig1tar := 1.4142 * as_float( split( private.mysens.calcsensitivity1('gain'))[1] ); snr1tar := targetflux[i] / sig1tar phase1tar := 4.06 * 10 / snr1tar; # this is the thermal noise --> phase in the estimate for the cross-band phase difference phase1all := sqrt( phase1tar^2 + (phase1cal*targetfreq/calfreq)^2 ); nints0[i] := as_integer( ( phase1all / gainerror )^2 + 0.5 ); if (nints0[i] <= 1) { nints0[i] := 1; } tgerr0[i] := phase1all / sqrt(nints0[i]); # now, get the total time used for this calibration tslew0[i] := private.slew( x[i]/cos(el/(180/pi) ), y[i] ); tint0[i] := nints0[i] * ( 2*tchangefreq + 2 ); toverhead0[i] := 2*tslew0[i] + tint0[i]; } # determine the optimal source and its consequences toverhead := toverhead0[( crosscycletime > toverhead0 )]; tslew := tslew0[( crosscycletime > toverhead0 )]; tint := tint0[( crosscycletime > toverhead0 )]; nints := nints0[( crosscycletime > toverhead0 )]; tgerr := tgerr0[ ( crosscycletime > toverhead0 )]; # there is also a small amount of phase error in the cross-band phase diff, # due to the 2 second cycle time atmophase := private.gamma0 * private.stidata.getphase( freq=targetfreq, baseline=v*2, el=el); tgerr2 := sqrt( tgerr^2 + atmophase^2 ); if (len(toverhead) == 0 ) { results.ok := F; } else { results.ok := T; # vectors of efficiencies for cross-band calibration eff1 := sqrt( (crosscycletime-toverhead)/crosscycletime ); eff2 := exp( - ((tgerr2/(180/pi))^2 / 2) ) eff0 := eff1 * eff2; iambest := order( - eff0 )[1]; results.eff0 := eff0[iambest]; results.eff1 := eff1[iambest]; results.duty := (crosscycletime-toverhead[iambest])/crosscycletime; results.eff2 := eff2[iambest]; results.optresid := tgerr2[iambest]; results.atmophase := atmophase; results.noisephase := tgerr[iambest] results.cbgainerror := results.atmophase / 3; results.calfluxtarget := targetflux[iambest]; results.calfluxcal := calflux[iambest]; results.caldistance := sqrt( x[iambest]^2 + y[iambest]^2) results.nints := nints[iambest]; results.toverhead := toverhead[iambest]; results.tslew := tslew[iambest]; results.tcal := tint[iambest]; } } return results; } # # based on dvt results & crossband results, identify the optimal cal strategy # self.evaluate := function( dvtrecord=[=], crossbandrecord=[=], targetfreq=345, el=60, gainerror=5 ) { evaldata := [=]; phaseresid := private.gamma0 * private.stidata.getphase( freq=targetfreq, baseline=dvtrecord.dvt, el=el) phaseresid1 := sqrt( phaseresid^2 + gainerror^2 ); # efficiency due to just time off source fs_eff1 := sqrt( (dvtrecord.beta -1 )/ dvtrecord.beta ); # efficiency due to just decorrelation fs_eff2 := e^( - ( (phaseresid1/(180/pi))^2 )/2 ); fs_eff0 := fs_eff1 * fs_eff2; iorder := order( -fs_eff0 ); evaldata.beta := dvtrecord.beta[iorder[1]]; # this will give the optimal beta evaldata.fs := [=]; evaldata.all := [=]; evaldata.fs.eff0 := fs_eff0[iorder[1]]; evaldata.fs.eff1 := fs_eff1[iorder[1]]; evaldata.fs.eff2 := fs_eff2[iorder[1]]; evaldata.fs.duty := (evaldata.beta -1 )/ evaldata.beta ; evaldata.all.eff0 := fs_eff0[iorder[1]] * crossbandrecord.eff0; evaldata.all.eff1 := fs_eff1[iorder[1]] * crossbandrecord.eff1; evaldata.all.eff2 := fs_eff2[iorder[1]] * crossbandrecord.eff2; evaldata.fs.optresid := phaseresid1[iorder[1]]; evaldata.fs.atmophase := phaseresid[iorder[1]]; evaldata.all.optresid := sqrt( evaldata.fs.optresid^2 + crossbandrecord.optresid^2); evaldata.fs.gainerror := gainerror; if (evaldata.fs.optresid > 4.5 * gainerror) { # gainerror is too tightly spec'ed for the atmos evaldata.fs.gainerror := evaldata.fs.optresid / 4.5; } else if (evaldata.fs.optresid < gainerror) { # gainerror is worse than the atmosphere! evaldata.fs.gainerror := evaldata.fs.optresid; } evaldata.fs.optdvt := dvtrecord.dvt[iorder[1]]; return evaldata; } # For a list of calibration sources, dvt estimates # - which calibrator is best # - integration time on calibrator # - optimal cycle time (considering phase decorrelation) # - optimum d + vt/2 # # strategy: calculate the d+vt/2 for tcycle = beta * toverhead # where beta := 0.5 * [4:40] # Another program with knowledge of the phase fluctuations can # then weigh the duty cycle noise with the decorrelation noise # # sourcelist [=] record of sources and distances # tchangefreq [s] time required to change frequencies # el [deg] elevation of observations # v [m/s] velocity of the atmosphere # gainerror [deg] specifie gain error (sets tcal) # beta1 dim-less initial guess for beta self.dvt := function(sourcelist0=[=], calfreq=90, targetfreq=345, tchangefreq=1.0, el=60, v=10, gainerror=5.0, beta1=8.0, maxdist=5) { # select out the inner part of the sourcelist sourcelist := [=]; x := sourcelist0.x; y := sourcelist0.y; sourcelist.freq := sourcelist0.freq; sourcelist.x := sourcelist0.x[ x > -maxdist & x < maxdist & y > -maxdist & y < maxdist ]; sourcelist.y := sourcelist0.y[ x > -maxdist & x < maxdist & y > -maxdist & y < maxdist ]; sourcelist.s := sourcelist0.s[ x > -maxdist & x < maxdist & y > -maxdist & y < maxdist ]; sourcelist.alpha := sourcelist0.alpha[ x > -maxdist & x < maxdist & y > -maxdist & y < maxdist ]; x := F; y := F; # OK, we need to consider things like the SNR we need on the calibrator # calibrator's gain error must be more accurate than the target # note: in radians now! gainerror := max( gainerror, 0.1); # guard against insanity # HEY! SNR of 10 ==> rms phase of 4.06 degrees.... but then we must # multiply by sqrt(2) as we are differencing two gains to get the # baseline-based phase == 5.74 degrees # snr_required := 10 * ( 5.74 / gainerror ) * targetfreq / calfreq; if (el >= 90) { el := 89; } airmass := 1/sin( el /(180/pi) ); # pick the calibrator we think should be best (we are guessing at beta) nsources := len(sourcelist.s); dvt0 := 999999; d0 := 999999; iopt := 1; toverhead0 := 999999; tcal0 := 999999; tslew0 := 999999; d0 := 999999; calflux := 0 * sourcelist.s; if (calfreq > 90) { calflux := sourcelist.s * (calfreq/sourcelist.freq)^ (-sourcelist.alpha - private.kinkalpha); } else { calflux := sourcelist.s * (calfreq/sourcelist.freq)^(-sourcelist.alpha); } for (i in [1:nsources]) { sigma := calflux[i] / snr_required; private.mysens.setatmos(tau225=private.tau225, airmass=airmass ); private.mysens.setobservation(npol=2, bandwidth=8, freq=calfreq, dt=1); tcal := private.mysens.calcinttime1(sigma=sigma, mode='gain'); tslew := private.slew( sourcelist.x[i]/cos(el/(180/pi) ), sourcelist.y[i] ); tmove := max(tslew, tchangefreq); toverhead := tcal + 2 * tmove; tcycle := toverhead * beta1; delta := sqrt( (sourcelist.x[i])^2 + (sourcelist.y[i])^2 ); d := delta/(180/pi) * private.h * airmass; dvt := (d + v * tcycle/2) ; if (dvt < dvt0) { d0 := d; dvt0 := dvt; toverhead0 := toverhead; tcal0 := tcal; tslew0 := tslew; iopt := i; } } myrec := [=]; myrec.dx := sourcelist.x[iopt]; myrec.dy := sourcelist.y[iopt]; myrec.s := calflux[iopt]; myrec.d := d0 myrec.tslew := tslew0; myrec.tcal := tcal0; myrec.toverhead := toverhead0; # Note: for different values of beta, different cals may be optimal myrec.beta := 0.5*[4:40]; tcycle := toverhead0 * myrec.beta; myrec.dvt := ( d + v * tcycle / 2) ; return ref myrec; } ########################################################### # # Private Functions # ########################################################### # slew - calculates the time it takes to slew in az/el # NOTE: you must supply the actual daz you want, pre-divided # by cos(el) # # model for fast switching: derived from VERTEX simulations # each move consists of a START, a SLEW, and a SETTLE # the time spent in each phase is parameterized in terms of the # distance, as fit to the 1.5 and 4 deg VERTEX simulations # # daz = required slew in azimuth [deg] # del = required slew in elevation [deg] # private.slew := function(daz=0, del=0) { daz := abs(daz); del := abs(del); tazstart := 0.152 + 0.052 * daz; tazstop := 0.430 + 0.080 * daz; azslope := min( (2 + .7 * daz), 6); tazslew := daz/azslope; taz := tazslew + tazstart + tazstop; telstart := 0.254 + 0.044 * del; telstop := 0.414 + 0.044 * del; elslope := min((1.724 + 0.224 * del), 3); telslew := del/elslope; tel := telslew + telstart + telstop; t := max(taz, tel); return t; } }