# this is where we estimate a template for the bandpass shape using # [calon(scan x) - calon(scan x+2) ] / calon(scan x+2) + # [caloff(scan x) - caloff(scan x+2) ] / caloff(scan x+2) # assume that the data have been calibrated and that the correct # Tsys values are in the headers baselinetemplate := function (val scana, val scanb, val nif = 1, val nfeed = 1) { if (nfeed == 1) { # get data for feed that goes on/off vsigon := d.getr(scana, phase=2, nfeed=1, nif=nif); # scan x, cal vsigoff := d.getr(scana, phase=1, nfeed=1, nif=nif); # scan x, nocal vrefon := d.getr(scanb, phase=2, nfeed=1, nif=nif); # scan x+2, cal vrefoff := d.getr(scanb, phase=1, nfeed=1, nif=nif); # scan x+2, nocal } else { vsigon := d.getr(scanb, phase=2, nfeed=2, nif=nif); # scan x+2, cal vsigoff := d.getr(scanb, phase=1, nfeed=2, nif=nif); # scan x+2, nocal vrefon := d.getr(scana, phase=2, nfeed=2, nif=nif); # scan x, cal vrefoff := d.getr(scana, phase=1, nfeed=2, nif=nif); # scan x, nocal } # get tsys values from the header - average tsys from the four scans tsys1r := (vsigon.header.tsys[1] + vsigoff.header.tsys[1] + vrefon.header.tsys[1] + vrefoff.header.tsys[1])/4; tsys1l := (vsigon.header.tsys[2] + vsigoff.header.tsys[2] + vrefon.header.tsys[2] + vrefoff.header.tsys[2])/4; # print " tsys=",tsys1r,tsys1l; # average together the cal on and cal off data vsig := vsigon; vsig.data.arr := (vsigoff.data.arr + vsigon.data.arr)/2.; vsig.data.arr := (vsigoff.data.arr + vsigon.data.arr)/2.; vref:= vrefon; vref.data.arr := (vrefoff.data.arr + vrefon.data.arr)/2.; # now calculate the [scan(x) - scan(x+2)]/scan(x+2) bandpass template vsig.data.arr[1,] := tsys1r*(vsig.data.arr[1,] - vref.data.arr[1,]) / vref.data.arr[1,]; vsig.data.arr[2,] := tsys1l*(vsig.data.arr[2,] - vref.data.arr[2,]) / vref.data.arr[2,]; # put the tsys value in the header vsig.header.tsys[1] := tsys1r; vsig.header.tsys[2] := tsys1l; # d.plotscan(vsig,T); # put the resulting data into the globalscan d.uniput("globalscan1",vsig); return T; } # calculate the bandpass template for a set of scans # scans is a list of scans containing the first scan of a nod baselinetemplateLoop := function (scans, val nif = 1) { d.sclear(); for (i in 1:(len(scans)-1)) { j := scans[i]; k := scans[i+1]; print "Processing scans", j,k," and",j+1,k+1; baselinetemplate(j, k, nif, 1); d.accum() baselinetemplate(j, k, nif, 2); d.accum() baselinetemplate(j+1, k+1, nif, 1); d.accum() baselinetemplate(j+1, k+1, nif, 2); d.accum() } d.ave(); d.show(); } # given the SD record of the data and the basline template calculate # a fit of data = a + b * frequency + c * baseline template # and then remove that fit and return resulting data # third input is the region (index) where the fit is to be made # fourth input is a mask of which values to fit # first the fitter must be included include 'functionfitter.g' baselinetemplatefit := function(sdData,sdTemplate,region,fitmask=[T,T,T]) { # will need oposite boolean values for the function fitter whattofit := !fitmask; # start the function fitter myfitter := functionfitter(); # data array size size := sdData.data.arr::shape; # resulting sdRecord to subtract from data fitbaseline := sdTemplate; # loop by polarization for (pol in 1:size[1]) { y := array(0,size[2]); x := array(0,size[2]*2); datamask := array(F,size[2]); datamask[region] := T; # get data into funky state needed by functionfitter y := sdData.data.arr[pol,]; for (j in 1:size[2]) { k := (j-1)*2+1; x[k] := j; x[k+1] := sdTemplate.data.arr[pol,j]; } # give the data to the function fitter myfitter.setdata(x,y,mask=datamask); # give the function to fit (x0=frequency x1=baseline template myfitter.setfunction('p0 + p1*x0 + p2*x1'); # set initial guess guess := [0.068, 0.0, 1.0]; myfitter.setparameters(guess); # now calculate the fit and get errors and chi square ans := myfitter.fit(linear = T, fixed=whattofit); anserr := myfitter.geterror(); anschisq := myfitter.getchisq(); print ans,anserr,anschisq; # now calculate baseline to remove fitbaseline.data.arr[pol,] *:= ans[3]; fitbaseline.data.arr[pol,] +:= ans[1]; fitbaseline.data.arr[pol,] +:= [1:size[2]]*ans[2]; 1.e9; } # subtract the fit from the data and plot sdData.data.arr -:= fitbaseline.data.arr; d.plotscan(sdData,F); # put into globalscan d.uniput("globalscan1",sdData); # close the fitter myfitter.done(); return T; } ## here are my commands # # d.open('kuACS') # scans := seq(25,155,2)[seq(25,155,2)!=75 & seq(25,155,2)!=31 & seq(25,155,2)!=67] # scans := seq(37,155,2)[seq(37,155,2)!=75 & seq(37,155,2)!=67] # ## calibrate data # for (s in scans) { # d.calib(s) # nod data now done correctly # } # # include '/users/tminter/scintle3/GBT_01A_004/nonlinear/baselinetemplate.g' # ## now average first spectral window # d.sclear() # for (s in scans) { # d.plotc(s,nif=1) # d.accum() # } # d.ave() # d.show() # intAvgSpect1 := d.gsget() # ## now find baseline template for first spectral window # baselinetemplateLoop(scans,1) # intBaseTemp1 := d.gsget() # ## plot data and baseline template # d.plotscan(intAvgSpect1,F) # d.plotscan(intBaseTemp1,T) # ## smooth the data ## d.plotscan(intAvgSpect1) ## d.boxcar(16) ## intAvgSpect1sm := d.gsget() ## d.plotscan(intBaseTemp1) ## d.boxcar(16) ## intBaseTemp1sm := d.gsget() # # fit the baseline template to the data # baselinetemplatefit(intAvgSpect1,intBaseTemp1,[125:1975]) # baselinetemplatefit(intAvgSpect1sm,intBaseTemp1sm,[1:32]) # # intCorrected1 := d.gsget() # ## now average first spectral window # d.sclear() # for (s in scans) { # d.plotc(s,nif=2) # d.accum() # } # d.ave() # d.show() # intAvgSpect2 := d.gsget() # ## now find baseline template for first spectral window # baselinetemplateLoop(scans,2) # intBaseTemp2 := d.gsget() # # plot data and baseline template # d.plotscan(intAvgSpect2,F) # d.plotscan(intBaseTemp2,T) # # fit the baseline template to the data # baselinetemplatefit(intAvgSpect2,intBaseTemp2,[125:1975]) # # intCorrected2 := d.gsget()