# Glish script to flag and calibrate a 334 MHz observation of W3
   # After calibration, the data are gridded into a spectral line cube.
   # HISTORY
   # 04Aug19 GIL compute S-R/R directly
   # 04Aug18 GIL test of calibrating using first and last scans
   # 04Aug17 GIL use marsTemp.g code
   # 04Aug17 GIL complete full calibration
   # 03Apr30 GIL update for m nord's run
   # 03Apr01 GIL Initial version
   #
   # This script assumes that the AIPS++ dish package has been loaded.
   # In GB this is done by typing:
   # % gbdish
   #             at the linux command prompt
   # utilities needed to do the calculations
   include '/users/glangsto/flareGlish/ott.g'
   include '/users/glangsto/flareGlish/polAvg.g'
   include 'reduce1B.g'

   offScan  := 705
   d.open('mapRef_sd')
   mapRef := d.getscan(offScan)
   d.plotscan(mapRef)
   refData := getData( mapRef)

   d.open('calJy3C286_sd')
   cal3c286JyK := d.getscan(512)
   d.plotscan( cal3c286JyK)
   calData := getData( cal3c286JyK)

   # Nord's observations on 2003 July 4th had two complete maps with 91 scans
   # RaLongMap: 523:613   ;  DecLatMap 614:704
   # Import scans into a measurement set for the AIPS++ dish package
   dataName :='/home/archive/science-data/tape-0002/AGBT02A_030_04'
   msName   :='AGBT02A_030_04_SP'
   firstScan:=614
   lastScan :=704
   d.import( dataName,,,firstScan,lastScan)
   # Examine the overall properties of the Measurement set
   d.gms()
   # define the numbers scans to examine and image
   scans := [firstScan:lastScan]

   print 'Tsys =', mapRef.header.tsys

   #for (i in scans) { d.SRcal( i, mapRef); d.plotc( i)}

   # scale calibration data for reference Tsys
   calData[1,] := calData[1,] * mapRef.header.tsys[1]
   calData[2,] := calData[2,] * mapRef.header.tsys[2]

   t := table(msName,readonly=F);
   allData := t.getcol('FLOAT_DATA')
   shape(allData)

   dataShape := shape( allData)
   nPol     := dataShape[1]
   nChan    := dataShape[2]
   nSpectra := dataShape[3]

   for (p in 1:nPol) {   # for both polarizations
     for (i in 1:nChan) { # for all channels in a spectrum
      allData[p,i,] := calData[p,i] * (allData[p,i,] - refData[p,i])/refData[p,i]
     }
   } # for all polarizations 

   t.putcol('CORRECTED_DATA',allData);
   t.done()

   # if flagging
   doFlag := F
   if (doFlag) {
    flags := [1:56,255:270,360:370,414:420,511:513,720:740,760:770,976:1024]
    for (i in scans) { 
      d.flag( i,,,flags,T);
    }
   } # end if flagging

   # if mapping....
   doMap := F
   if (doMap) { 
    include 'image.g'
    imageName := 'sgrL.im'
    d.imagems(msName,100,900,10,imname=imageName,cellx='6arcmin',celly='6arcmin',nx=250,ny=250,gridfn='PB')
    sgrA := image(imageName)
   #set the beam size in arc seconds
    sgrA.setrestoringbeam(major=40.4*60,minor=40.4*60,pa=0)
    sgrA.setbrightnessunit('Jy/beam')
   }