# # pragma include once; include 'image.g' include 'quanta.g' include 'measures.g' include 'pgplotter.g' include 'otfpath.g' include 'fitting.g' include 'almasensitivity.g' include 'phasemonitor.g' include 'randomnumbers.g'; include 'oneoverfnoise.g'; exploreelevation := function() { rec30:= optimalvslew(freq=300, sourcesize=.2, pwv=1.0, stiphase=1.0, el=30 ); rec40:= optimalvslew(freq=300, sourcesize=.2, pwv=1.0, stiphase=1.0, el=40 ); rec50:= optimalvslew(freq=300, sourcesize=.2, pwv=1.0, stiphase=1.0, el=50 ); rec60:= optimalvslew(freq=300, sourcesize=.2, pwv=1.0, stiphase=1.0, el=60 ); rec70:= optimalvslew(freq=300, sourcesize=.2, pwv=1.0, stiphase=1.0, el=70 ); print rec40; print rec50; print rec60; print rec70; return T; } # # stand-alone function which calculates the slew velocity which results # in the lowest sigma (including both residual atmosphere and thermal noise) # optimalvslew := function(freq=300, sourcesize=.2, pwv=1.0, stiphase=1.0, el=60, az=0.0, fitorder=3, gainerr=0.0 ) { nscans:=15; mysim := sdatmosim(); mysim.readaips('ATMO.DOUBLE') vslews := [ .1, .2, .4, .6, .8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0]; # # the FAST list # vslews := [.1, .4, .8, 1.2]; sigma := 0 * vslews; dtvec := 0 * vslews; tscan := 0 * vslews; dutycycle := 0 * vslews; #print 'DEBUG!!! in optimalvslew: freq = ', freq ; for (ii in ind(vslews)) { vslew := vslews[ii]; dt := mysim.suggestdt(freq=freq, vslew=vslew, el=el) print 'Suggested dt = ', dt; mysim.setup(elevation=el, azimuth=az, pwv=pwv, stiphase=stiphase, freq=freq, gainfluctuations=gainerr); mysim.setupotf1(sourcesize=sourcesize, dt=dt, nscans=nscans, vslew=vslew); if (ii == 1) { mysim.viewsetup(); } results1 := mysim.otf1(addnoise=T, fitorder=fitorder); print results1; sigma[ii] := results1.rms * sqrt( results1.tscan ); dtvec[ii] := dt; tscan[ii] := results1.tscan; dutycycle[ii] := results1.dutycycle; } results := [=]; for (ii in ind(vslews)) { print vslews[ii], sigma[ii] } results.vslew := vslews; results.sigma := sigma; results.tscan := tscan; results.dutycycle := dutycycle; results.bestvslew := sort_pair( sigma, vslews )[1] results.bestsigma := sort( sigma )[1] results.freq := freq; results.el := el; results.az := az; results.pwv := pwv; results.stiphase := stiphase; results.sourcesize := sourcesize; # now, compare with the atmosphere-less case: dt := mysim.suggestdt(freq=freq, vslew=results.bestvslew, el=el) mysim.setup(azimuth=90, elevation=el, pwv=pwv, stiphase=0.000001, freq=freq, gainfluctuations=gainerr); mysim.setupotf1(sourcesize=sourcesize, dt=dt, nscans=nscans, vslew=results.bestvslew); results1 := mysim.otf1(addnoise=T, fitorder=2); results.bestdt := dt; results.bestsigmanoatm := results1.rms * sqrt( results1.tscan ); return results; } # # define the key tool for doing these simulations # sdatmosim := subsequence() { private := [=]; private.atmosaips := ''; # name of AIPS++ atmo image private.atmosname1:= ''; # name of FITS atmo image private.atmos1 := F; # current atmospheric image - don't use private.arr := F; # current atmospheric array - use this private.imdelta := 0.5; # m private.shape := []; private.dishdiam := 12.0; # m private.velocity := 12.0; # m/s private.height := 500.0; # m private.temperature:= 250.0; # Kelvins private.pwv := 1.0; # mm private.stiphase := 1.0; # [deg] phase @ 300m, 11.198 GHz private.alpha := F; # wet tau = alpha * pwv private.freq := 300.0; # GHz private.starttime := 0.0; # starting time, seconds private.az := F; # azimuth, rad private.el := F; # el, rad private.bs1data := [=]; # holding data for bs1 private.bs1data.dt:= F; # integration tims of bs1 private.bs1data.tset:= F; # setup time for bs1 private.bs2data := [=]; # holding data for bs2 private.bs2data.dt:= F; # integration tims of bs2 private.bs2data.tset:= F; # setup time for bs2 private.bs3data := [=]; # holding data for bs3 private.bs3data.dt:= F; # integration tims of bs3 private.bs3data.tset:= F; # setup time for bs3 private.otf1data := [=]; # holding data for otf private.otf1data.dt:= F; # integration tims of otf1 private.sti.baseline := 300; # site testing interferometer bl [m] private.sti.freq := 11.198e+9; # sti freq [Hz] private.sti.c := 3e+8; # [m/s] private.sti.factor := 6.3/1000; # path [m] caused by 1mm PWV private.pg := F; # pg plotter tool private.mysens := almasensitivity(); # sensitivity tool private.mysens.setgaincor(F); # turn off e^{tau.airmass} correction private.phmon := phasemonitor(); # phase monitor tool private.phmon.setsite('alma'); private.gainfluc := 0.0001; # 1 part in 10^4 1/f noise private.seed := 1005; # random number seed private.oneoverf := F; # 1/f noise generator tool private.tsysbackground := F; # tsys background level, required for 1/f # and reset by resetbackground(); # public functions self.type := function() { # note ('I do single dish atmospheric noise simulations'); return 'sdatmosim'; } self.resetbackground := function() { wider private; private.tsysbackground := F; return T; } self.getprivate := function() { return ref private; } self.done := function() { wider self, private; if (is_tool(private.atmos1)) private.atmos1.done(); if (is_tool(private.mysens)) private.mysens.done(); if (is_tool(private.phmon)) private.phmon.done(); if (is_tool(private.oneoverf)) private.oneoverf.done(); private.donepgplotter(); self := F; val public := F; return T; } # Set up paramaters and coordinates for the Beam Switching Simulation # az, el in degrees; throw in arcmin; dt in seconds == HALF CYCLE TIME # tsetup in seconds: how much time is lost moving # dt is the total cycle time self.setupbs1 := function(mythrow=1, dt=0.05, tset=0.005, ncycles=200) { wider private; if (!private.checkatmos()) { print 'Before you can run setupbs1(), you'; print 'need to get an ATMOS image with readaips() or readfits()'; return F; } if (private.el == F) { print 'setupbs1() requires setup() first'; return F; } if (tset >= dt) { print 'HEY! dt is less than tset! fix that in setupbs1()'; return F; } private.bs1data := [=]; private.bs1data.throw := mythrow/ 60 / 57.2957795; # now in rad private.bs1data.dt := dt; private.bs1data.tset := tset; private.bs1data.ncycles := ncycles; private.bs1data.t := private.starttime + dt * [1:ncycles]; private.beampixels := private.dishdiam/private.imdelta/2.0 # setup coordinates OFF ON OFF ON doffset := (private.bs1data.throw * private.height) / private.imdelta / sin( private.el ); xoffset := doffset * cos( private.az ); yoffset := doffset * sin( private.az ); print 'xoffset, yoffset = ', xoffset, yoffset; private.bs1data.ix1 := private.shape[1] - as_integer(1.5*private.beampixels) - xoffset; private.bs1data.iy1 := as_integer(private.shape[2]/2); private.bs1data.xon := private.bs1data.ix1 - (private.bs1data.t + dt/2) * private.velocity/private.imdelta; private.bs1data.yon := private.bs1data.iy1 + 0 * private.bs1data.t; private.bs1data.on := 0 * private.bs1data.t; private.bs1data.xoff := private.bs1data.ix1 + xoffset - private.bs1data.t * private.velocity/private.imdelta; private.bs1data.yoff := private.bs1data.iy1 + yoffset + 0 * private.bs1data.t; private.bs1data.off := 0 * private.bs1data.t; # determine if we are going to run into trouble ixlast := private.bs1data.ix1 - ncycles * dt * private.velocity/private.imdelta; if (ixlast > (1.5*private.beampixels)) { print 'Sims will be fine spatially'; return T; } else { print 'ATMO Model is too small to accomodate this long a sim'; print 'ixlast = ', ixlast; return F; } } # OFF1-ON1-OFF1-ON1-OFF1-ON1 OFF2-ON1-OFF2-ON1-OFF2-ON1 # nsmallcycles=3 # \--------------------- nbigcycles = 2 ------------------/ # total OFFS (or ONS): nsmallcycles * 2 * nbigcycles # Set up paramaters and coordinates for Traditional Double Beam Switching self.setupbs2 := function(mythrow=1, dt=0.05, tset=0.005, nsmallcycles=10, nbigcycles=20) { wider private; # nbigcycles must be even nbigcycles := 2 * as_integer((nbigcycles+1)/2); nsmallcycles := as_integer(nsmallcycles); if (!private.checkatmos()) { print 'Before you can run setupbs2(), you'; print 'need to get an ATMOS image with readaips() or readfits()'; return F; } if (private.el == F) { print 'setupbs1() requires setup() first'; return F; } if (tset >= dt) { print 'HEY! dt is less than tset! fix that in setupbs2()'; return F; } private.beampixels := private.dishdiam/private.imdelta/2.0; throwrad := mythrow / 60 / 57.2957795; doffset := (throwrad * private.height) / private.imdelta / sin( private.el); xoffset := doffset * cos( private.az ); yoffset := doffset * sin( private.az); tstart := private.starttime; # # Note: because of the way we access the data here (private.bs2data[i]), # each of these indexed records must be defined BEFORE any other data members # of private.bs2data (such as private.bs2data.nbigcycles) # private.bs2data := [=]; for (i in [1:nbigcycles]) { private.bs2data[i] := [=]; private.bs2data[i].t := tstart + dt * [1:nsmallcycles]; private.bs2data[i].dt := dt; tstart +:= dt * nsmallcycles; if (i == 2*as_integer(i/2.0)) { isign := 1; } else { isign := -1; } private.bs2data[i].ix1 := private.shape[1] - as_integer(1.5*private.beampixels) - xoffset; private.bs2data[i].iy1 := as_integer(private.shape[2]/2); # note: 0 * t ==> makes it a vector of the right length private.bs2data[i].xon := private.bs2data[i].ix1 - (private.bs2data[i].t + dt/2) * private.velocity/private.imdelta; private.bs2data[i].yon := private.bs2data[i].iy1 + 0 * private.bs2data[i].t; private.bs2data[i].on := 0 * private.bs2data[i].t; private.bs2data[i].xoff := private.bs2data[i].ix1 + isign*xoffset - private.bs2data[i].t * private.velocity/private.imdelta; private.bs2data[i].yoff := private.bs2data[i].iy1 + isign*yoffset + 0 * private.bs2data[i].t; private.bs2data[i].off := 0 * private.bs2data[i].t; } # # finished defining indexed members of private.bs2data # private.bs2data.nbigcycles := nbigcycles; private.bs2data.nsmallcycles := nsmallcycles; private.bs2data.throw := throwrad; private.bs2data.dt := dt; private.bs2data.tset := tset; private.beampixels := private.dishdiam/private.imdelta/2.0 # determine if we are going to run into trouble ixlast := private.bs2data[1].ix1 - nsmallcycles *nbigcycles * dt * private.velocity/private.imdelta; if (ixlast > (1.5*private.beampixels)) { print 'BS2 Sims will be fine spatially'; return T; } else { print 'BS2: ATMO Model is too small to accomodate this long a sim'; print 'ixlast = ', ixlast; return F; } } # Set up paramaters and coordinates for the Beam Switching Simulation # az, el in degrees; throw in arcmin; dt in seconds # dt is the total cycle time self.setupbs3 := function(mythrow=1, dt=0.05, tset=0.005, ncycles=200) { wider private; if (!private.checkatmos()) { print 'Before you can run setupbs3(), you'; print 'need to get an ATMOS image with readaips() or readfits()'; return F; } if (private.el == F) { print 'setupbs3() requires setup() first'; return F; } if (tset >= dt) { print 'HEY! dt is less than tset! fix that in setupbs3()'; return F; } private.bs3data := [=]; private.bs3data.throw := mythrow/ 60 / 57.2957795; # now in rad private.bs3data.dt := dt; private.bs3data.tset := tset; private.bs3data.ncycles := ncycles; private.bs3data.t := private.starttime + dt * [1:ncycles]; private.beampixels := private.dishdiam/private.imdelta/2.0 # setup coordinates OFF1 ON OFF2 ON doffset := (private.bs3data.throw * private.height) / private.imdelta / sin( private.el ); xoffset := doffset * cos( private.az ); yoffset := doffset * sin( private.az ); print 'xoffset, yoffset = ', xoffset, yoffset; private.bs3data.ix1 := private.shape[1] - as_integer(1.5*private.beampixels) - xoffset; private.bs3data.iy1 := as_integer(private.shape[2]/2); private.bs3data.xon := private.bs3data.ix1 - (private.bs3data.t + dt/2) * private.velocity/private.imdelta; private.bs3data.yon := private.bs3data.iy1 + 0 * private.bs3data.t; private.bs3data.on := 0 * private.bs3data.t; isign := [1:len(private.bs3data.t)]; jsign := 0 * isign; jsign[(isign % 2) == 0] := -1; jsign[(isign % 2) == 1] := 1; private.bs3data.xoff := private.bs3data.ix1 + jsign * xoffset - private.bs3data.t * private.velocity/private.imdelta; private.bs3data.yoff := private.bs3data.iy1 + jsign * yoffset + 0 * private.bs3data.t; private.bs3data.off := 0 * private.bs3data.t; # determine if we are going to run into trouble ixlast := private.bs3data.ix1 - ncycles * dt * private.velocity/private.imdelta; if (ixlast > (1.5*private.beampixels)) { print 'Sims will be fine spatially'; return T; } else { print 'ATMO Model is too small to accomodate this long a sim'; print 'ixlast = ', ixlast; return F; } } # now run the BS1 beam switching simulation # return results as a record self.bs1 := function(interp=T, addnoise=F) { wider private; results := [=]; if (is_boolean(private.bs1data.dt)) { print 'Please run setupbs1() before running bs1()'; return F; } results := private.bs0( private.bs1data, interp, addnoise ); results.observation := 'beam_switch_1'; results.dutycycle := 0.5; results.dt := private.bs1data.dt; results.throw := private.bs1data.throw; return ref results; } # This one may not do anything different! # now run the BS2 beam switching simulation # return results as a record self.bs2 := function(interp=T, addnoise=F) { wider private; results := [=]; if (is_boolean(private.bs2data.dt)) { print 'Please run setupbs2() before running bs2()'; return F; } radiuspix := ( private.dishdiam / 2.0 )/ private.imdelta; diff := []; nn := 0; results.rmsresid := 0*[1:(private.bs2data.nbigcycles/2)]; for (i in [1:private.bs2data.nbigcycles]) { oneresult := private.bs0( private.bs2data[i], interp, addnoise, returndiff=T ); # OK: we do this elaborate thing here to cancel out atmospheric trends; # This is like traditional double beam switching # The sqrt(2) factor in the rms calculation is to correct for the # depletion of random noise by adding the two diff vector if (i != (2*as_integer(i/2)) ) { # 1st of each pair: diff := oneresult.diff; } else { # second of each pair diff +:= oneresult.diff; results.rmsresid[i/2] := sqrt( 2 * sum(diff^2)/len(diff) ); } } results.rms := sum( results.rmsresid ) / len( results.rmsresid ); results.dutycycle := 0.5; results.observation := 'beam_switch_2'; if (interp) { results.method := 'interpolated'; } else { results.method := 'difference'; } results.dt := private.bs2data.dt; results.throw := private.bs2data.throw; return ref results; } # now run the BS3 beam switching simulation # return results as a record self.bs3 := function(interp=T, addnoise=F) { wider private; results := [=]; if (is_boolean(private.bs3data.dt)) { print 'Please run setupbs3() before running bs3()'; return F; } results := private.bs0( private.bs3data, interp, addnoise ); results.observation := 'beam_switch_3'; results.dutycycle := 0.5; results.dt := private.bs3data.dt; results.throw := private.bs3data.throw; return ref results; } # sourcesize in deg # velocity in deg/sec # dt in seconds self.setupotf1 := function(sourcesize=0.1, vslew=0.5, jerkmax=100, dt=0.016, nscans=10) { wider private; if (!private.checkatmos()) { print 'Before you can run setupotf1(), you'; print 'need to get an ATMOS image with readaips() or readfits()'; return F; } if (private.el == F) { print 'setupotf1() requires setup() first'; return F; } private.otf1data := F; # # Note: because of the way we access the data here (private.otf1data[i]), # each of these indexed records must be defined BEFORE any other data members # of private.otf1data (such as private.otf1data.nscans) # mypath := otfpath(); # OK: for a source 1 deg, at el=60deg, we need to scan 2 deg in az mypath.setup(v=vslew, x=sourcesize/cos(private.el), jmax=jerkmax, dt=0.001); # info on time-on-source and time to turn around are contant in 'results' # also note: by docalc(*T*), we include the full turnaround in the trajectory, # but not in the duty cycle and all. results := mypath.docalc(T); tscan := results.tonsource + results.toffsource; private.beampixels := private.dishdiam/private.imdelta/2.0 iy1 := as_integer(private.shape[2]/2); ix1 := private.shape[1] - as_integer(1.5*private.beampixels) - private.height * sourcesize/57.2957795/sin(private.el); private.otf1data := [=]; toffset := -(results.toffsource)/2; # will this fit or not? trajectory := mypath.gettrajectory(dt=dt, forward=T); islew := (private.height/sin(private.el)) * (trajectory.x*cos(private.el)/57.295775) / private.imdelta; ix := ix1 - (trajectory.t + toffset) * private.velocity/ private.imdelta + cos(private.az) * islew; dix := (ix[len(ix)] - ix[1]) * (len(ix)/(len(ix)-1)); dixn := dix * nscans; xfinal := ix[1] + dixn; print 'DIX, DIXN, XFINAL: ', dix, dixn, xfinal if (xfinal < 20) { print 'This is trouble, we are running off the edge of the ATMOS' nscans := as_integer( abs((ix[1] -20 ) / dix )); print 'So, we are changing nscans to', nscans } forward := ( ([1:nscans] % 2) != 0 ); for (i in [1:nscans]) { private.otf1data[i] := [=]; trajectory := mypath.gettrajectory(dt=dt, forward=forward[i]); private.otf1data[i].t := trajectory.t + toffset; toffset := private.otf1data[i].t[len(private.otf1data[i].t)] - (results.toffsource)/2; # listen carefully: height/sin(el) == distance to the turbulent layer # trajectory.x * cos(el) accounts for the fact that # the antenna needs to move 2 degrees to cover 1 deg # at el=60deg; /57.29 turns it into radians , mult # by hight/sin(el) converts the angle to an arc in the # turbulent layer # /imdelta turns it into a pixel number islew := ( private.height / sin(private.el) ) * (trajectory.x * cos(private.el) /57.295775) / private.imdelta; private.otf1data[i].ix := ix1 - private.otf1data[i].t * private.velocity/private.imdelta + cos(private.az) * islew private.otf1data[i].iy := iy1 + sin(private.az) * islew private.otf1data[i].onoroff := trajectory.phase; private.otf1data[i].sky := 0.0 * trajectory.phase; } private.otf1data.nscans := nscans; private.otf1data.sourcesize := sourcesize; private.otf1data.vslew := vslew; private.otf1data.jerkmax := jerkmax; private.otf1data.dutycycle := results.dutycycle; private.otf1data.dt := dt; private.otf1data.tscan := tscan; private.otf1data.nave := self.beamfwhm(private.freq) / (2 * 60 * vslew * cos( private.el) * dt ); print 'Scan Details: tscan = ', tscan; print ' : duty = ', results.dutycycle; print ' : length= ', len(private.otf1data[1].onoroff); print ' : nscans= ', nscans; print ' Averaging WOULD be: ', private.otf1data.nave; return T; } # this el is in degrees! self.suggestdt := function(freq=300, vslew=0.5, el=60 ) { #print 'DEBUG!!! in suggestdt: freq = ', freq ; dt := self.beamfwhm(freq) / (2 * 60 * vslew * cos(el/57.3)); return dt; } # perform an OTF observation # We will have a region to one side called OFF1, a region on the # other called OFF2, and the ON measurements; we # perform a poly fit to the OFF1 and OFF2 regions and interpolate # to remove from the ON values # Hey, we also need to average within each beam (small effect?) self.otf1 := function(addnoise=F, fitorder=1) { wider private; if (is_boolean(private.otf1data.dt)) { print 'Please run setupotf1() before running otf1()'; return F; } fitorder := as_integer(fitorder); if (fitorder < 1) { fitorder := 1; print 'Setting fitorder to 1 to get average atmosphere out'; } results := [=]; results.rmsresid := 0 * [1:private.otf1data.nscans]; radiuspix := ( private.dishdiam / 2.0 )/ private.imdelta; noiserms := 0.0; rand:=randomnumbers(); rand.reseed(); if (addnoise ) { private.mysens.setdt(private.otf1data.dt); noiserms2 := private.mysens.calcsensitivity1('obs'); note(paste('We will add thermal noise of', noiserms2)); noiserms := dq.quantity(noiserms2).value; } for (i in [1:private.otf1data.nscans]) { for (j in ind(private.otf1data[i].t)) { private.otf1data[i].sky[j] := private.observe( private.otf1data[i].ix[j], private.otf1data[i].iy[j], radiuspix); } # add thermal noise if (addnoise) { noise := rand.normal(mean=0.0, variance=(noiserms^2), shape=len(private.otf1data[i].sky)); private.otf1data[i].sky +:= noise; } # multiply by 1/f-infected gains if (private.gainfluc > 0) { private.applygain (private.otf1data[i].t, private.otf1data.dt, private.otf1data[i].sky) } # fit baseline to atmosphere myfit := fitter(); xx := private.otf1data[i].t[(private.otf1data[i].onoroff == 0)]; xxave := sum(xx)/len(xx); xx := xx - xxave; yy := private.otf1data[i].sky[(private.otf1data[i].onoroff == 0)]; myfit.init(fitorder); myfit.makepoly(xx, yy); myfit.fit(); sol := myfit.solution(); ont := private.otf1data[i].t[(private.otf1data[i].onoroff == 1)] - xxave; onsky:= private.otf1data[i].sky[(private.otf1data[i].onoroff == 1)]; fitsky := 0 * onsky; for (k in [1:fitorder]) { fitsky +:= sol[k] * ont^(k-1) } onsky -:= fitsky; results.rmsresid[i] := sqrt( sum ( onsky^2 ) / (len(onsky)-1) / sin(private.el) ); # airmass effect) } rand.done(); # scale up by e^{tau*airmass}, the gain correction tau := private.mysens.getwet(private.freq) * private.pwv + private.mysens.getdry(private.freq); results.rmsresid *:= exp( tau * (1/sin(private.el)) ); results.rms := sum(results.rmsresid)/len(results.rmsresid); #was 1/sqrt( sum(1/(results.rmsresid)^2 ) ) results.observation := 'on_the_fly_1'; results.fitorder := fitorder; results.noise := noiserms; results.tscan := private.otf1data.tscan; results.nscans := private.otf1data.nscans; results.sourcesize := private.otf1data.sourcesize; results.vslew := private.otf1data.vslew; results.jerkmax := private.otf1data.jerkmax; results.dutycycle := private.otf1data.dutycycle; results.dt := private.otf1data.dt; print 'OTF data sims'; return results; } self.gettsys := function() { if (private.el == F) { print 'gettsys() requires setup() first'; return F; } tsys := private.mysens.calctsys(); return tsys; } self.beamfwzi := function(freq=300) { lambda := 0.3 / freq; beamrads := 2 * 1.22 * lambda / private.dishdiam; beamminutes := beamrads * 180 / pi * 60; return beamminutes; } self.beamfwhm := function(freq=300) { lambda := 0.3 / freq; beamrads := 1.22 * lambda / private.dishdiam; beamminutes := beamrads * 180 / pi * 60; return beamminutes; } self.summary := function() { return (self.viewsetup()); } self.viewsetup := function() { print 'AZ [deg] = ', private.az * 57.2957795; print 'EL [deg] = ', private.el * 57.2957795; print 'V [m/s] = ', private.velocity; print 'Diam [m] = ', private.dishdiam; print 'H [m] = ', private.height; print 'temp [K] = ', private.temperature; print 'freq [GHz]= ', private.freq; print '1/f flucs = ', private.gainfluc; note(paste('AZ [deg] = ', private.az * 57.2957795)); note(paste('EL [deg] = ', private.el * 57.2957795)); note(paste('V [m/s] = ', private.velocity)); note(paste('Diam [m] = ', private.dishdiam)); note(paste('H [m] = ', private.height)); note(paste('temp [K] = ', private.temperature)); note(paste('freq [GHz]= ', private.freq)); note(paste('alpha = ', private.alpha)); note(paste('1/f flucs = ', private.gainfluc)); return T; } self.setup := function(azimuth=0, elevation=60, diam=12, velocity=12, height=500, tatmos=250, freq=300, pwv=1, stiphase=1.0, gainfluctuations=0.0001, seed = 1005) { wider private; if (elevation <= 0 || elevation > 90) { print 'Please use an elevation > 0 deg and <= 90 deg'; return F; } private.pwv := pwv; if (stiphase < 1e-5) { stiphase := 1e-5; } private.stiphase := stiphase; private.az := azimuth/57.2957795; private.el := elevation/57.2957795; private.dishdiam := diam; private.velocity := abs(velocity); if (velocity < 0) { print 'Our geometrical assumptions require a positive velocity'; print 'It has been made positive for your convenience'; } private.height := height; private.temperature := tatmos; private.freq := freq; private.starttime := 0.0; private.mysens.setobservation(npol=2, bandwidth=4, freq=freq, dt=1); private.mysens.setinstrument (nants=1, diam=diam, baseline=diam); private.mysens.setatmos (tatmos=tatmos, tamb=(tatmos+20), pwv=pwv, airmass=(1/sin(private.el)) ); private.mysens.usekelvins(); private.alpha := private.mysens.getwet(private.freq); if (is_tool(private.atmos1)) { # needed to update private.arr for PWV and STIPHASE private.mungearr(); } private.gainfluc := gainfluctuations; private.seed := seed; if (is_tool(private.oneoverf)) private.oneoverf.done(); if (private.gainfluc > 0) { private.oneoverf := oneoverfnoise(); # 10368 is composite and will not easily repeat for runovers # 10368 * 0.005 = 83 sec # # for what its worth: we found the flattest Allan Variances with a # power law of 1.05 # private.oneoverf.setup(n=10368, dt=0.008, power=1.05, seed=private.seed); private.oneoverf.makenoise(targetmean=1.0, targetrms=private.gainfluc, timerms=1); } self.resetbackground(); return T; } self.readaips := function(atmoaips='') { wider private; if (is_tool(private.atmos1)) { private.atmos1.done(); } private.atmosaips :=atmoaips; private.atmos1 := image(private.atmosaips); private.getimageinfo(); if (!is_fail(private.atmos1)) { return T; } else { print 'Could not read file ', atmoaips; return F; } } self.readfits := function(atmosfits='', istrip=1, blc=[1,1], trc=[2048,300]) { wider self, private; private.atmosaips := spaste( atmosfits, '.aips') private.atmosname1 := spaste( atmosfits, '.strip', istrip) im0 := imagefromfits(infile=atmosfits, outfile=private.atmosaips); myshape := im0.shape(); trc[1] := min(myshape[1], trc[1]); trc[2] := min(myshape[2], trc[2]); blc[1] := max(blc[1], 1); blc[2] := max(blc[2], 1); blc1 := myshape; trc1 := myshape; trc1[1] := trc[1]; trc1[2] := trc[2]; blc1[1] := blc[1]; blc1[2] := blc[2]; r1 := drm.box(blc1, trc1); private.atmos1 := im0.subimage(region=r1, outfile=private.atmosname1); im0.done(); private.getimageinfo(); return 'maybe'; } self.view := function() { if (private.checkatmos()) { print 'Note that this image does not reflect PWV rescaling'; private.atmos1.view(); return T; } else { return F; } } # returns the structure function as a record: # results.xx[] # results.rms[] # results.fitamp # results.fitexp # results.avepwv # results.stiphase self.getspatialstructurefunc := function(doplot=T) { return (private.phmon.getspatialstructurefunc(doplot)); } #==================== private functions ================================================ # get coords info, atmos array, scale pixel values for PWV private.getimageinfo := function() { wider private; if (!private.checkatmos()) { print 'No image found in getimageinfo()'; fail; } cs := private.atmos1.coordsys(); private.imdelta := cs.i()[1]; cs.done(); private.arr := private.atmos1.getchunk(); private.shape := private.atmos1.shape(); private.mungearr(); return T; } # set the PWV AVE and RMS to agree with what we want private.mungearr := function() { ### now munge the arr so that it gives the correct PWV and # STIPHASE ################################################ # make sure phmon tool has the array! private.phmon.setarray( private.arr, private.imdelta ); structurefunction := self.getspatialstructurefunc(doplot=F); # make it zero mean; rms unchanged private.arr -:= structurefunction.avepwv; # scale fluctuations to desired level private.arr *:= private.stiphase / structurefunction.stiphase; # still zero mean: now add a constant level to get PWV true private.arr +:= private.pwv; private.phmon.setarray( private.arr, private.imdelta ); ### munge complete ####################################### return T; } # Beam Switching primitive, reused by several BS sims private.bs0 := function( ref base=[=], interp=T, addnoise=F, returndiff=F) { radiuspix := ( private.dishdiam / 2.0 )/ private.imdelta; if (len(base.t) <= 1) { print 'base.t not found in bs0() primitive'; fail; } noiserms := 0.0; rand:=randomnumbers(); rand.reseed(); # note that base.dt is for the ON-OFF cycle time; # for the noise calc, we use base.dt - base.tset if (addnoise) { private.mysens.setdt(base.dt - base.tset); noiserms2 := private.mysens.calcsensitivity1('obs'); note(paste('We will add thermal noise of', noiserms2)); noiserms := dq.quantity(noiserms2).value; } results := [=]; for (i in ind(base.t)) { base.on[i] := private.observe(base.xon[i], base.yon[i], radiuspix); base.off[i] := private.observe(base.xoff[i], base.yoff[i], radiuspix); } if (addnoise) { noise := rand.normal(mean=0.0, variance=(noiserms^2), shape=len(base.on)); base.on +:= noise; noise := rand.normal(mean=0.0, variance=(noiserms^2), shape=len(base.off)); base.off +:= noise; } # multiply by 1/f-infected gains if (private.gainfluc > 0) { private.applygain (base.t, base.dt/2, base.off) private.applygain ((base.t + base.dt/2), base.dt/2, base.on) } if (interp) { # print 'Using interpolation to determine beam_switching_1 residual'; nsamples := len(base.on); ON := base.on[1:(nsamples-1)]; OFF1 := base.off[1:(nsamples-1)]; OFF2 := base.off[2:nsamples]; diff := ON - (OFF1 + OFF2)/2; ON := F; OFF1 := F; OFF2 := F; results.method := 'interpolated'; } else { # print 'No interpolation to determine beam_switching_1 residual'; diff := base.on - base.off results.method := 'difference'; } tau := private.mysens.getwet(private.freq) * private.pwv + private.mysens.getdry(private.freq); results.rms := sqrt( sum(diff^2)/ len(diff) / sin(private.el) ) * exp( tau * (1/sin(private.el)) ); if (returndiff) { results.diff := diff * exp( tau * (1/sin(private.el)) ); } diff := F return ref results; } # get the average gain over each integration and multiply the # data elements by the gains private.applygain := function(timevec=[], dt=0.01, ref datavec=[]) { for (i in ind(timevec)) { gain := private.oneoverf.average(t=timevec[i], dt=dt); datavec[i] *:= gain; } return T; } 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.checkatmos := function() { if (is_fail(private.atmos1) || is_boolean(private.atmos1) || private.atmos1.type() != 'image') { print 'No image found in getimageinfo()'; return F; } else { return T; } } # Given the atmosphere, the beam center, and the beam column, # calculate the beam atmospheric contribution. # This calculation DOES NOT include elevation effects! # # Hey, we added tsysbackground for use with 1/f calcs # private.observe := function(xcenter, ycenter, radiuspix) { wider private; if (!private.checkatmos()) { print 'In observe(), problem with private.atmos1'; return F; } blc := private.shape; trc := private.shape; blc[1] := as_integer(xcenter) - radiuspix; trc[1] := as_integer(xcenter) + radiuspix + 1; blc[2] := as_integer(ycenter) - radiuspix; trc[2] := as_integer(ycenter) + radiuspix + 1; if (blc[1] < 1) { print 'observe(): blc x problems'; fail; } if (blc[2] < 1) { print 'observe(): blc y problems'; fail; } if (trc[1] > private.shape[1]) { print 'observe(): trc x problems'; fail; } if (trc[2] > private.shape[2]) { print 'observe(): trc y problems'; fail; } xnewcenter := xcenter - as_integer(xcenter) + radiuspix + 1; ynewcenter := ycenter - as_integer(ycenter) + radiuspix + 1; atmosarr := private.arr[ blc[1]:trc[1], blc[2]:trc[2],,,]; beamarr := atmosarr * 0 + 1 dist2arr := atmosarr * 0 rad2 := radiuspix^2; ix := [1:(2*radiuspix + 2)] dx2 := (ix - xnewcenter)^2 for (iy in [1:(2*radiuspix + 2)]) { dist2arr[ix, iy] := dx2 + (iy - ynewcenter)^2; } beamarr[(dist2arr > rad2)] := 0.0; ntotal := sum(beamarr); skytotal := sum( beamarr * atmosarr ); atmosarr := F; beamarr := F; dist2arr := F; if (ntotal > 0) { result := private.temperature * (1-exp(-private.alpha * skytotal/ntotal)); if (is_boolean(private.tsysbackground) ) { private.tsysbackground := private.mysens.calctsys() - result; } return (result + private.tsysbackground); } else { print 'No pixels selected in observe()'; print 'xcenter=',xcenter,'ycenter=',ycenter; print 'xnewcenter=',xnewcenter,'ynewcenter=',ynewcenter; print 'rad2', rad2, 'blc', blc, 'trc', trc; fail; } } } # end of constructor compare2atmos := function() { include 'sdatmosim.g' mysim := sdatmosim() print 'Old Atmos' mysim.readaips('ATMO.DOUBLE') mysim.setup(azimuth=180, elevation=60, pwv=2, stiphase=5.0, freq=300) mysim.getspatialstructurefunc(); mysim.setupbs1(mythrow=0, dt=0.1, ncycles=100) results1 := mysim.bs1(F) results1 results1 := mysim.bs1(T) results1 print 'New Atmos' mysim.readaips('ATMO.DOUBLE') mysim.setup(azimuth=180, elevation=60, pwv=2, stiphase=5.0, freq=300) mysim.getspatialstructurefunc(); mysim.setupbs1(mythrow=0, dt=0.1, ncycles=100) mysim.bs1(F) mysim.bs1(T) mysim.setupbs1(mythrow=0, dt=0.2, ncycles=60) mysim.bs1(F) mysim.bs1(T) } simpletest := function() { include 'sdatmosim.g' mysim := sdatmosim() mysim.readaips('ATMO.DOUBLE') mysim.getspatialstructurefunc(); freq := 300; el := 60; vslew := 0.5 mysim.setup(azimuth=225, elevation=el, pwv=1, stiphase=2.0, freq=freq, gainfluctuations=0.0) dt := mysim.suggestdt(freq=freq, vslew=vslew, el=el) mysim.setupotf1(sourcesize=0.1, vslew=vslew, jerkmax=100, dt=dt, nscans=10) results1 := mysim.otf1(addnoise=F, fitorder=3); results2 := mysim.otf1(addnoise=T, fitorder=3); print results1; print results2; print 'Tsys = ',(mysim.gettsys()), 'K'; mysim.done(); return T; } gaintest := function() { include 'sdatmosim.g' mysim := sdatmosim() mysim.readaips('ATMO.DOUBLE') mysim.getspatialstructurefunc(); freq := 300; el := 60; vslew := 0.5 # no NOTHING mysim.setup(azimuth=225, elevation=el, pwv=1, stiphase=0.0000001, freq=freq, gainfluctuations=0.0) dt := mysim.suggestdt(freq=freq, vslew=vslew, el=el) mysim.setupotf1(sourcesize=0.06, vslew=vslew, jerkmax=100, dt=dt, nscans=1) results1 := mysim.otf1(addnoise=F, fitorder=1); # 3.32151561e-08, results2 := mysim.otf1(addnoise=F, fitorder=2); # 1.94249816e-08 results3 := mysim.otf1(addnoise=F, fitorder=3); # 4.49566914e-09 # just THERMAL NOISE results1 := mysim.otf1(addnoise=T, fitorder=1); # 0.0171899752 results2 := mysim.otf1(addnoise=T, fitorder=2); # 0.0158172198 results3 := mysim.otf1(addnoise=T, fitorder=3); # 0.0168764766 # 10-3 gains mysim.setup(azimuth=225, elevation=el, pwv=1, stiphase=0.0000001, freq=freq, gainfluctuations=0.001) dt := mysim.suggestdt(freq=freq, vslew=vslew, el=el) mysim.setupotf1(sourcesize=0.06, vslew=vslew, jerkmax=100, dt=dt, nscans=1) results1 := mysim.otf1(addnoise=F, fitorder=1); results2 := mysim.otf1(addnoise=F, fitorder=2); results3 := mysim.otf1(addnoise=F, fitorder=3); # 10-3 gains plus noise results1 := mysim.otf1(addnoise=T, fitorder=1); results2 := mysim.otf1(addnoise=T, fitorder=2); results3 := mysim.otf1(addnoise=T, fitorder=3); print 'Tsys = ',(mysim.gettsys()), 'K'; mysim.done(); return T; } testsdatmosim := function() { include 'sdatmosim.g' mysim := sdatmosim() mysim.readaips('ATMO.DOUBLE') mysim.getspatialstructurefunc(); mysim.setup(azimuth=225, elevation=60, pwv=1, stiphase=2.0, freq=300) mysim.getspatialstructurefunc(); mysim.setupbs1(mythrow=1, dt=0.05, ncycles=200) results1 := mysim.bs1(T) results1 mysim.setup(azimuth=225, elevation=60, pwv=1, freq=220) results1 := mysim.bs1(T) results1 mysim.setup(azimuth=225, elevation=60, pwv=1, freq=90) results1 := mysim.bs1(T) results1 mysim.setupbs2(mythrow=1, dt=0.05, nsmallcycles=10, nbigcycles=10) results2 := mysim.bs2(F) mysim.setup(azimuth=225, elevation=60, pwv=1, freq=300); mysim.setupbs3(mythrow=1, dt=0.05, ncycles=200) results1 := mysim.bs3(T) mysim.setup(azimuth=225, elevation=60, pwv=1, freq=90) mysim.setupotf1(sourcesize=0.2, vslew=0.5, jerkmax=100, dt=0.016, nscans=10) mysim.otf1(addnoise=F, fitorder=2); mysim.otf1(addnoise=F, fitorder=4); mysim.otf1(addnoise=F, fitorder=6); mysim.otf1(addnoise=T, fitorder=2); mysim.otf1(addnoise=T, fitorder=4); mysim.otf1(addnoise=T, fitorder=6); mysim.setupotf1(sourcesize=0.2, vslew=0.5, jerkmax=100, dt=0.016, nscans=10) mysim.setup(azimuth=225, elevation=60, pwv=1, freq=345) mysim.otf1(addnoise=F, fitorder=2); mysim.otf1(addnoise=F, fitorder=4); mysim.otf1(addnoise=F, fitorder=6); mysim.otf1(addnoise=T, fitorder=2); mysim.otf1(addnoise=T, fitorder=4); mysim.otf1(addnoise=T, fitorder=6); mysim.done(); }