Build: #146 was successful
Job: Build and Package RHEL6 was successful
Build log
The build generated 9,234 lines of output.The output is too long and has been truncated to the last 1,000 lines. Download or view full build log.
| 18-Sep-2019 23:19:13 | + [pwvtime, pwv, antenna] = readpwv('%s/'%vis) |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("No CalWVR.xml file found, so no PWV retrieved. Copy it to this directory and try again.") |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + return(pwvmean,-1) |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + matches = np.where(np.array(pwvtime)>myTimes[0])[0] |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("Found no times > %d" % (myTimes[0])) |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + return(0,-1) |
| 18-Sep-2019 23:19:13 | + if (len(pwv) < 1): |
| 18-Sep-2019 23:19:13 | + print("Found no PWV data") |
| 18-Sep-2019 23:19:13 | + return(0,-1) |
| 18-Sep-2019 23:19:13 | + ptime = np.array(pwvtime)[matches] |
| 18-Sep-2019 23:19:13 | + matchedpwv = np.array(pwv)[matches] |
| 18-Sep-2019 23:19:13 | + matches2 = np.where(ptime<myTimes[-1])[0] |
| 18-Sep-2019 23:19:13 | + if (len(matches2) < 1): |
| 18-Sep-2019 23:19:13 | + # look for the value with the closest start time |
| 18-Sep-2019 23:19:13 | + mindiff = 1e12 |
| 18-Sep-2019 23:19:13 | + for i in range(len(pwvtime)): |
| 18-Sep-2019 23:19:13 | + if (abs(myTimes[0]-pwvtime[i]) < mindiff): |
| 18-Sep-2019 23:19:13 | + mindiff = abs(myTimes[0]-pwvtime[i]) |
| 18-Sep-2019 23:19:13 | + pwvmean = pwv[i]*1000 |
| 18-Sep-2019 23:19:13 | + matchedpwv = [] |
| 18-Sep-2019 23:19:13 | + for i in range(len(pwvtime)): |
| 18-Sep-2019 23:19:13 | + if (abs(abs(myTimes[0]-pwvtime[i]) - mindiff) < 1.0): |
| 18-Sep-2019 23:19:13 | + matchedpwv.append(pwv[i]) |
| 18-Sep-2019 23:19:13 | + pwvmean = 1000*np.median(matchedpwv) |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Taking the median of %d pwv measurements from all antennas = %.3f mm" % (len(matchedpwv),pwvmean)) |
| 18-Sep-2019 23:19:13 | + pwvstd = np.std(matchedpwv) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + pwvmean = 1000*np.median(matchedpwv[matches2]) |
| 18-Sep-2019 23:19:13 | + pwvstd = np.std(matchedpwv[matches2]) |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Taking the median of %d pwv measurements from all antennas = %.3f mm" % (len(matches2),pwvmean)) |
| 18-Sep-2019 23:19:13 | +# mytb.done() |
| 18-Sep-2019 23:19:13 | + return(pwvmean,pwvstd) |
| 18-Sep-2019 23:19:13 | +# end of getMedianPWV |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def computeAzElFromRADecMJD(raDec, mjd, observatory='ALMA'): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Computes the az/el for a specified J2000 RA/Dec, MJD and observatory. |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + raDec must be in radians: [ra,dec] |
| 18-Sep-2019 23:19:13 | + mjd must be in days |
| 18-Sep-2019 23:19:13 | + returns the [az,el] in radians |
| 18-Sep-2019 23:19:13 | + - Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + myme = createCasaTool(metool) |
| 18-Sep-2019 23:19:13 | + myqa = createCasaTool(qatool) |
| 18-Sep-2019 23:19:13 | + mydir = myme.direction('J2000', myqa.quantity(raDec[0],'rad'), myqa.quantity(raDec[1],'rad')) |
| 18-Sep-2019 23:19:13 | + myme.doframe(myme.epoch('mjd', myqa.quantity(mjd, 'd'))) |
| 18-Sep-2019 23:19:13 | + myme.doframe(myme.observatory(observatory)) |
| 18-Sep-2019 23:19:13 | + myazel = myme.measure(mydir,'azel') |
| 18-Sep-2019 23:19:13 | + myqa.done() |
| 18-Sep-2019 23:19:13 | + myme.done() |
| 18-Sep-2019 23:19:13 | + return([myazel['m0']['value'], myazel['m1']['value']]) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def getRADecForField(ms, myfieldId, debug): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Returns RA,Dec in radians for the specified field in the specified ms. |
| 18-Sep-2019 23:19:13 | + -- Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + myms = createCasaTool(mstool) |
| 18-Sep-2019 23:19:13 | + myms.open(ms) |
| 18-Sep-2019 23:19:13 | + myd = myms.getfielddirmeas('DELAY_DIR', fieldid=myfieldId) # dircolname defaults to 'PHASE_DIR' |
| 18-Sep-2019 23:19:13 | + myms.close() |
| 18-Sep-2019 23:19:13 | + mydir = np.array([[myd['m0']['value']], [myd['m1']['value']]]) # simulates tb.getcell |
| 18-Sep-2019 23:19:13 | + return(mydir) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def findClosestTime(mytimes, mytime): |
| 18-Sep-2019 23:19:13 | + myindex = 0 |
| 18-Sep-2019 23:19:13 | + mysep = np.abs(mytimes[0]-mytime) |
| 18-Sep-2019 23:19:13 | + for m in range(1,len(mytimes)): |
| 18-Sep-2019 23:19:13 | + if (np.abs(mytimes[m] - mytime) < mysep): |
| 18-Sep-2019 23:19:13 | + mysep = np.abs(mytimes[m] - mytime) |
| 18-Sep-2019 23:19:13 | + myindex = m |
| 18-Sep-2019 23:19:13 | + return(myindex) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def getWeather(vis='', scan='', antenna='0',verbose=False, mymsmd=None): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Queries the WEATHER and ANTENNA tables of an .ms by scan number or |
| 18-Sep-2019 23:19:13 | + list of scan numbers in order to return mean values of: angleToSun, |
| 18-Sep-2019 23:19:13 | + pressure, temperature, humidity, dew point, wind speed, wind direction, |
| 18-Sep-2019 23:19:13 | + azimuth, elevation, solarangle, solarelev, solarazim. |
| 18-Sep-2019 23:19:13 | + If the sun is below the horizon, the solarangle returned is negated. |
| 18-Sep-2019 23:19:13 | + -- Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Entered getWeather with vis,scan,antenna = %s,%s,%s" % (str(vis), str(scan), str(antenna))) |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + if str(antenna).isdigit(): |
| 18-Sep-2019 23:19:13 | + antennaName = mymsmd.antennanames(antenna)[0] |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + antennaName = antenna |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + antenna = mymsmd.antennaids(antennaName)[0] |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + antennaName = string.upper(antenna) |
| 18-Sep-2019 23:19:13 | + antenna = mymsmd.antennaids(antennaName)[0] |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("Either the ANTENNA table does not exist or antenna %s does not exist" % (antenna)) |
| 18-Sep-2019 23:19:13 | + return([0,[]]) |
| 18-Sep-2019 23:19:13 | + mytb = createCasaTool(tbtool) |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + mytb.open("%s/POINTING" % vis) |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("POINTING table does not exist") |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + return([0,0]) |
| 18-Sep-2019 23:19:13 | + subtable = mytb.query("ANTENNA_ID == %s" % antenna) |
| 18-Sep-2019 23:19:13 | + mytb.close() |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + mytb.open("%s/OBSERVATION" % vis) |
| 18-Sep-2019 23:19:13 | + observatory = mytb.getcell("TELESCOPE_NAME",0) |
| 18-Sep-2019 23:19:13 | + mytb.close() |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("OBSERVATION table does not exist, assuming observatory == ALMA") |
| 18-Sep-2019 23:19:13 | + observatory = "ALMA" |
| 18-Sep-2019 23:19:13 | + if (scan == ''): |
| 18-Sep-2019 23:19:13 | + scan = mymsmd.scannumbers() |
| 18-Sep-2019 23:19:13 | + conditions = {} |
| 18-Sep-2019 23:19:13 | + conditions['pressure']=conditions['temperature']=conditions['humidity']=conditions['dewpoint']=conditions['windspeed']=conditions['winddirection'] = 0 |
| 18-Sep-2019 23:19:13 | + conditions['scan'] = scan |
| 18-Sep-2019 23:19:13 | + if (type(scan) == str): |
| 18-Sep-2019 23:19:13 | + if (scan.find('~')>0): |
| 18-Sep-2019 23:19:13 | + tokens = scan.split('~') |
| 18-Sep-2019 23:19:13 | + scan = [int(k) for k in range(int(tokens[0]),int(tokens[1])+1)] |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + scan = [int(k) for k in scan.split(',')] |
| 18-Sep-2019 23:19:13 | + if (type(scan) == type(np.ndarray(0))): |
| 18-Sep-2019 23:19:13 | + scan = list(scan) |
| 18-Sep-2019 23:19:13 | + if (type(scan) == list): |
| 18-Sep-2019 23:19:13 | + myTimes = np.array([]) |
| 18-Sep-2019 23:19:13 | + for sc in scan: |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + print("calling timesforscan") |
| 18-Sep-2019 23:19:13 | + newTimes = mymsmd.timesforscan(sc) |
| 18-Sep-2019 23:19:13 | + print("times = %s" % (str(newTimes))) |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("Error reading scan %d, is it in the data?" % (sc)) |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + return([conditions,[]]) |
| 18-Sep-2019 23:19:13 | + myTimes = np.concatenate((myTimes,newTimes)) |
| 18-Sep-2019 23:19:13 | + elif (scan != None): |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + myTimes = mymsmd.timesforscan(scan) |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("Error reading scan %d, is it in the data?" % (scan)) |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + return([conditions,[]]) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + return([conditions,[]]) |
| 18-Sep-2019 23:19:13 | + if (type(scan) == str): |
| 18-Sep-2019 23:19:13 | + scan = [int(k) for k in scan.split(',')] |
| 18-Sep-2019 23:19:13 | + if (type(scan) == list): |
| 18-Sep-2019 23:19:13 | + listscan = "" |
| 18-Sep-2019 23:19:13 | + listfield = [] |
| 18-Sep-2019 23:19:13 | + for sc in scan: |
| 18-Sep-2019 23:19:13 | +# print("Processing scan ", sc) |
| 18-Sep-2019 23:19:13 | + listfield.append(mymsmd.fieldsforscan(sc)) |
| 18-Sep-2019 23:19:13 | + listscan += "%d" % sc |
| 18-Sep-2019 23:19:13 | + if (sc != scan[-1]): |
| 18-Sep-2019 23:19:13 | + listscan += "," |
| 18-Sep-2019 23:19:13 | +# print("listfield = ", listfield) |
| 18-Sep-2019 23:19:13 | + listfields = np.unique(listfield[0]) |
| 18-Sep-2019 23:19:13 | + listfield = "" |
| 18-Sep-2019 23:19:13 | + for field in listfields: |
| 18-Sep-2019 23:19:13 | + listfield += "%s" % field |
| 18-Sep-2019 23:19:13 | + if (field != listfields[-1]): |
| 18-Sep-2019 23:19:13 | + listfield += "," |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + listscan = str(scan) |
| 18-Sep-2019 23:19:13 | + listfield = mymsmd.fieldsforscan(scan) |
| 18-Sep-2019 23:19:13 | + [az,el] = ComputeSolarAzElForObservatory(myTimes[0], mymsmd) |
| 18-Sep-2019 23:19:13 | + [az2,el2] = ComputeSolarAzElForObservatory(myTimes[-1], mymsmd) |
| 18-Sep-2019 23:19:13 | + azsun = np.mean([az,az2]) |
| 18-Sep-2019 23:19:13 | + elsun = np.mean([el,el2]) |
| 18-Sep-2019 23:19:13 | + direction = subtable.getcol("DIRECTION") |
| 18-Sep-2019 23:19:13 | + azeltime = subtable.getcol("TIME") |
| 18-Sep-2019 23:19:13 | + subtable.close() |
| 18-Sep-2019 23:19:13 | + telescopeName = mymsmd.observatorynames()[0] |
| 18-Sep-2019 23:19:13 | + if (len(direction) > 0 and telescopeName.find('VLA') < 0 and telescopeName.find('NRO') < 0): |
| 18-Sep-2019 23:19:13 | + azimuth = direction[0][0]*180.0/math.pi # a list of values |
| 18-Sep-2019 23:19:13 | + elevation = direction[1][0]*180.0/math.pi # a list of values |
| 18-Sep-2019 23:19:13 | + npat = np.array(azeltime) |
| 18-Sep-2019 23:19:13 | + matches = np.where(npat>myTimes[0])[0] |
| 18-Sep-2019 23:19:13 | + matches2 = np.where(npat<myTimes[-1])[0] |
| 18-Sep-2019 23:19:13 | + if (len(matches2) > 0 and len(matches) > 0): |
| 18-Sep-2019 23:19:13 | + if verbose: print("matches[0]=%d, matches2[-1]=%d" % (matches[0],matches[-1])) |
| 18-Sep-2019 23:19:13 | + matchingIndices = range(matches[0],matches2[-1]+1) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + matchingIndices = [] |
| 18-Sep-2019 23:19:13 | + if (len(matchingIndices) > 0): # CAS-8440 |
| 18-Sep-2019 23:19:13 | + conditions['azimuth'] = np.mean(azimuth[matches[0]:matches2[-1]+1]) |
| 18-Sep-2019 23:19:13 | + conditions['elevation'] = np.mean(elevation[matches[0]:matches2[-1]+1]) |
| 18-Sep-2019 23:19:13 | + elif (len(matches) > 0): # CAS-8440 |
| 18-Sep-2019 23:19:13 | + if verbose: print("using mean of all az/el values after time 0") |
| 18-Sep-2019 23:19:13 | + conditions['azimuth'] = np.mean(azimuth[matches[0]]) |
| 18-Sep-2019 23:19:13 | + conditions['elevation'] = np.mean(elevation[matches[0]]) |
| 18-Sep-2019 23:19:13 | + else: # CAS-8440 |
| 18-Sep-2019 23:19:13 | + if verbose: print("using mean of all az/el values") |
| 18-Sep-2019 23:19:13 | + conditions['azimuth'] = np.mean(azimuth) |
| 18-Sep-2019 23:19:13 | + conditions['elevation'] = np.mean(elevation) |
| 18-Sep-2019 23:19:13 | + conditions['solarangle'] = angularSeparation(azsun,elsun,conditions['azimuth'],conditions['elevation']) |
| 18-Sep-2019 23:19:13 | + conditions['solarelev'] = elsun |
| 18-Sep-2019 23:19:13 | + conditions['solarazim'] = azsun |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Using antenna = %s to retrieve mean azimuth and elevation" % (antennaName)) |
| 18-Sep-2019 23:19:13 | + print("Separation from sun = %f deg" % (abs(conditions['solarangle']))) |
| 18-Sep-2019 23:19:13 | + if (elsun<0): |
| 18-Sep-2019 23:19:13 | + conditions['solarangle'] = -conditions['solarangle'] |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Sun is below horizon (elev=%.1f deg)" % (elsun)) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Sun is above horizon (elev=%.1f deg)" % (elsun)) |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Average azimuth = %.2f, elevation = %.2f degrees" % (conditions['azimuth'],conditions['elevation'])) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + if (verbose): print("The POINTING table is either blank or does not contain Azim/Elev.") |
| 18-Sep-2019 23:19:13 | + if (type(scan) == int or type(scan)==np.int32): |
| 18-Sep-2019 23:19:13 | + # compute Az/El for this scan |
| 18-Sep-2019 23:19:13 | + myfieldId = mymsmd.fieldsforscan(scan) |
| 18-Sep-2019 23:19:13 | + if (type(myfieldId) == list or type(myfieldId) == type(np.ndarray(0))): |
| 18-Sep-2019 23:19:13 | + myfieldId = myfieldId[0] |
| 18-Sep-2019 23:19:13 | + fieldName = mymsmd.namesforfields(myfieldId) |
| 18-Sep-2019 23:19:13 | + if (type(fieldName) == list or type(fieldName) == type(np.ndarray(0))): |
| 18-Sep-2019 23:19:13 | + fieldName = fieldName[0] |
| 18-Sep-2019 23:19:13 | +# print("A) fieldname = ", fieldName) |
| 18-Sep-2019 23:19:13 | +# print("myfieldId = ", myfieldId) |
| 18-Sep-2019 23:19:13 | + myscantime = np.mean(mymsmd.timesforscan(scan)) |
| 18-Sep-2019 23:19:13 | +# print("Calling getRADecForField") |
| 18-Sep-2019 23:19:13 | + mydirection = getRADecForField(vis, myfieldId, verbose) |
| 18-Sep-2019 23:19:13 | + if (verbose): print("mydirection= %s" % (str(mydirection))) |
| 18-Sep-2019 23:19:13 | + if (len(telescopeName) < 1): |
| 18-Sep-2019 23:19:13 | + telescopeName = 'ALMA' |
| 18-Sep-2019 23:19:13 | + myazel = computeAzElFromRADecMJD(mydirection, myscantime/86400., telescopeName) |
| 18-Sep-2019 23:19:13 | + conditions['elevation'] = myazel[1] * 180/math.pi |
| 18-Sep-2019 23:19:13 | + conditions['azimuth'] = myazel[0] * 180/math.pi |
| 18-Sep-2019 23:19:13 | + conditions['solarangle'] = angularSeparation(azsun,elsun,conditions['azimuth'],conditions['elevation']) |
| 18-Sep-2019 23:19:13 | + conditions['solarelev'] = elsun |
| 18-Sep-2019 23:19:13 | + conditions['solarazim'] = azsun |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Separation from sun = %f deg" % (abs(conditions['solarangle']))) |
| 18-Sep-2019 23:19:13 | + if (elsun<0): |
| 18-Sep-2019 23:19:13 | + conditions['solarangle'] = -conditions['solarangle'] |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Sun is below horizon (elev=%.1f deg)" % (elsun)) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Sun is above horizon (elev=%.1f deg)" % (elsun)) |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Average azimuth = %.2f, elevation = %.2f degrees" % (conditions['azimuth'],conditions['elevation'])) |
| 18-Sep-2019 23:19:13 | + elif (type(scan) == list): |
| 18-Sep-2019 23:19:13 | + myaz = [] |
| 18-Sep-2019 23:19:13 | + myel = [] |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Scans to loop over = %s" % (str(scan))) |
| 18-Sep-2019 23:19:13 | + for s in scan: |
| 18-Sep-2019 23:19:13 | + fieldName = mymsmd.fieldsforscan(s) |
| 18-Sep-2019 23:19:13 | + if (type(fieldName) == list): |
| 18-Sep-2019 23:19:13 | + # take only the first pointing in the mosaic |
| 18-Sep-2019 23:19:13 | + fieldName = fieldName[0] |
| 18-Sep-2019 23:19:13 | + myfieldId = mymsmd.fieldsforname(fieldName) |
| 18-Sep-2019 23:19:13 | + if (type(myfieldId) == list or type(myfieldId)==type(np.ndarray(0))): |
| 18-Sep-2019 23:19:13 | + # If the same field name has two IDs (this happens in EVLA data) |
| 18-Sep-2019 23:19:13 | + myfieldId = myfieldId[0] |
| 18-Sep-2019 23:19:13 | + myscantime = np.mean(mymsmd.timesforscan(s)) |
| 18-Sep-2019 23:19:13 | + mydirection = getRADecForField(vis, myfieldId, verbose) |
| 18-Sep-2019 23:19:13 | + telescopeName = mymsmd.observatorynames()[0] |
| 18-Sep-2019 23:19:13 | + if (len(telescopeName) < 1): |
| 18-Sep-2019 23:19:13 | + telescopeName = 'ALMA' |
| 18-Sep-2019 23:19:13 | + myazel = computeAzElFromRADecMJD(mydirection, myscantime/86400., telescopeName) |
| 18-Sep-2019 23:19:13 | + myaz.append(myazel[0]*180/math.pi) |
| 18-Sep-2019 23:19:13 | + myel.append(myazel[1]*180/math.pi) |
| 18-Sep-2019 23:19:13 | + conditions['azimuth'] = np.mean(myaz) |
| 18-Sep-2019 23:19:13 | + conditions['elevation'] = np.mean(myel) |
| 18-Sep-2019 23:19:13 | + conditions['solarangle'] = angularSeparation(azsun,elsun,conditions['azimuth'],conditions['elevation']) |
| 18-Sep-2019 23:19:13 | + conditions['solarelev'] = elsun |
| 18-Sep-2019 23:19:13 | + conditions['solarazim'] = azsun |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Using antenna = %s to retrieve mean azimuth and elevation" % (antennaName)) |
| 18-Sep-2019 23:19:13 | + print("Separation from sun = %f deg" % (abs(conditions['solarangle']))) |
| 18-Sep-2019 23:19:13 | + if (elsun<0): |
| 18-Sep-2019 23:19:13 | + conditions['solarangle'] = -conditions['solarangle'] |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Sun is below horizon (elev=%.1f deg)" % (elsun)) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Sun is above horizon (elev=%.1f deg)" % (elsun)) |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Average azimuth = %.2f, elevation = %.2f degrees" % (conditions['azimuth'],conditions['elevation'])) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + # now, get the weather |
| 18-Sep-2019 23:19:13 | + if not os.path.exists('%s/WEATHER' % vis): |
| 18-Sep-2019 23:19:13 | + print("There is no WEATHER table for this ms.") |
| 18-Sep-2019 23:19:13 | + if (needToClose_mymsmd): mymsmd.close() |
| 18-Sep-2019 23:19:13 | + return([conditions,myTimes]) |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + mytb.open("%s/WEATHER" % vis) |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("Could not open the WEATHER table for this ms.") |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + return([conditions,myTimes]) |
| 18-Sep-2019 23:19:13 | + if (True): |
| 18-Sep-2019 23:19:13 | + mjdsec = mytb.getcol('TIME') |
| 18-Sep-2019 23:19:13 | + indices = np.argsort(mjdsec) |
| 18-Sep-2019 23:19:13 | + mjd = mjdsec/86400. |
| 18-Sep-2019 23:19:13 | + pressure = mytb.getcol('PRESSURE') |
| 18-Sep-2019 23:19:13 | + relativeHumidity = mytb.getcol('REL_HUMIDITY') |
| 18-Sep-2019 23:19:13 | + temperature = mytb.getcol('TEMPERATURE') |
| 18-Sep-2019 23:19:13 | + if (np.mean(temperature) > 100): |
| 18-Sep-2019 23:19:13 | + # must be in units of Kelvin, so convert to C |
| 18-Sep-2019 23:19:13 | + temperature -= 273.15 |
| 18-Sep-2019 23:19:13 | + if 'DEW_POINT' in mytb.colnames(): |
| 18-Sep-2019 23:19:13 | + dewPoint = mytb.getcol('DEW_POINT') |
| 18-Sep-2019 23:19:13 | + if (np.mean(dewPoint) > 100): |
| 18-Sep-2019 23:19:13 | + # must be in units of Kelvin, so convert to C |
| 18-Sep-2019 23:19:13 | + dewPoint -= 273.15 |
| 18-Sep-2019 23:19:13 | + if (np.mean(dewPoint) == 0): |
| 18-Sep-2019 23:19:13 | + # assume it is not measured and use NOAA formula to compute from humidity: |
| 18-Sep-2019 23:19:13 | + dewPoint = ComputeDewPointCFromRHAndTempC(relativeHumidity, temperature) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + dewPoint = None # Nobeyama measurement sets do not have a dewpoint column |
| 18-Sep-2019 23:19:13 | + sinWindDirection = np.sin(mytb.getcol('WIND_DIRECTION')) |
| 18-Sep-2019 23:19:13 | + cosWindDirection = np.cos(mytb.getcol('WIND_DIRECTION')) |
| 18-Sep-2019 23:19:13 | + windSpeed = mytb.getcol('WIND_SPEED') |
| 18-Sep-2019 23:19:13 | + mytb.done() |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + # put values into time order (they mostly are, but there can be small differences) |
| 18-Sep-2019 23:19:13 | + mjdsec = np.array(mjdsec)[indices] |
| 18-Sep-2019 23:19:13 | + pressure = np.array(pressure)[indices] |
| 18-Sep-2019 23:19:13 | + relativeHumidity = np.array(relativeHumidity)[indices] |
| 18-Sep-2019 23:19:13 | + temperature = np.array(temperature)[indices] |
| 18-Sep-2019 23:19:13 | + if dewPoint is not None: |
| 18-Sep-2019 23:19:13 | + dewPoint = np.array(dewPoint)[indices] |
| 18-Sep-2019 23:19:13 | + windSpeed = np.array(windSpeed)[indices] |
| 18-Sep-2019 23:19:13 | + sinWindDirection = np.array(sinWindDirection)[indices] |
| 18-Sep-2019 23:19:13 | + cosWindDirection = np.array(cosWindDirection)[indices] |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + # find the overlap of weather measurement times and scan times |
| 18-Sep-2019 23:19:13 | + matches = np.where(mjdsec>=np.min(myTimes))[0] |
| 18-Sep-2019 23:19:13 | + matches2 = np.where(mjdsec<=np.max(myTimes))[0] |
| 18-Sep-2019 23:19:13 | +# print("len(matches)=%d, len(matches2)=%d" % (len(matches), len(matches2))) |
| 18-Sep-2019 23:19:13 | + noWeatherData = False |
| 18-Sep-2019 23:19:13 | + if (len(matches)>0 and len(matches2) > 0): |
| 18-Sep-2019 23:19:13 | + # average the weather points enclosed by the scan time range |
| 18-Sep-2019 23:19:13 | + selectedValues = range(matches[0], matches2[-1]+1) |
| 18-Sep-2019 23:19:13 | + if (selectedValues == []): |
| 18-Sep-2019 23:19:13 | + # there was a either gap in the weather data, or an incredibly short scan duration |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("---- Finding the nearest weather value --------------------------- ") |
| 18-Sep-2019 23:19:13 | + selectedValues = findClosestTime(mjdsec, myTimes[0]) |
| 18-Sep-2019 23:19:13 | + elif (len(matches)>0): |
| 18-Sep-2019 23:19:13 | + # all points are greater than myTime, so take the first one |
| 18-Sep-2019 23:19:13 | + selectedValues = matches[0] |
| 18-Sep-2019 23:19:13 | + elif (len(matches2)>0): |
| 18-Sep-2019 23:19:13 | + # all points are less than myTime, so take the last one |
| 18-Sep-2019 23:19:13 | + selectedValues = matches2[-1] |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + # table has no weather data! |
| 18-Sep-2019 23:19:13 | + noWeatherData = True |
| 18-Sep-2019 23:19:13 | + if (noWeatherData): |
| 18-Sep-2019 23:19:13 | + conditions['pressure'] = 563.0 |
| 18-Sep-2019 23:19:13 | + conditions['temperature'] = 0 # Celsius is expected |
| 18-Sep-2019 23:19:13 | + conditions['humidity'] = 20.0 |
| 18-Sep-2019 23:19:13 | + conditions['dewpoint'] = -20.0 |
| 18-Sep-2019 23:19:13 | + conditions['windspeed'] = 0 |
| 18-Sep-2019 23:19:13 | + conditions['winddirection'] = 0 |
| 18-Sep-2019 23:19:13 | + print("WARNING: No weather data found in the WEATHER table!") |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + if (type(selectedValues) == np.int64 or type(selectedValues) == np.int32 or |
| 18-Sep-2019 23:19:13 | + type(selectedValues) == np.int): |
| 18-Sep-2019 23:19:13 | + conditions['readings'] = 1 |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("selectedValues=%d, myTimes[0]=%.0f, myTimes[1]=%.0f, len(matches)=%d, len(matches2)=%d" % (selectedValues, |
| 18-Sep-2019 23:19:13 | + myTimes[0],myTimes[1], len(matches), len(matches2))) |
| 18-Sep-2019 23:19:13 | + if (len(matches) > 0): |
| 18-Sep-2019 23:19:13 | + print("matches[0]=%f, matches[-1]=%f" % (matches[0], matches[-1])) |
| 18-Sep-2019 23:19:13 | + if (len(matches2) > 0): |
| 18-Sep-2019 23:19:13 | + print("matches2[0]=%f, matches2[-1]=%d" % (matches2[0], matches2[-1])) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + conditions['readings'] = len(selectedValues) |
| 18-Sep-2019 23:19:13 | + conditions['pressure'] = np.mean(pressure[selectedValues]) |
| 18-Sep-2019 23:19:13 | + if (conditions['pressure'] != conditions['pressure']): |
| 18-Sep-2019 23:19:13 | + # A nan value got through, due to no selected values (should be impossible)" |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print(">>>>>>>>>>>>>>>>>>>>>>>> selectedValues = %s" % (str(selectedValues))) |
| 18-Sep-2019 23:19:13 | + print("len(matches)=%d, len(matches2)=%d" % (len(matches), len(matches2))) |
| 18-Sep-2019 23:19:13 | + print("matches[0]=%f, matches[-1]=%f, matches2[0]=%f, matches2[-1]=%d" % (matches[0], matches[-1], matches2[0], matches2[-1])) |
| 18-Sep-2019 23:19:13 | + conditions['temperature'] = np.mean(temperature[selectedValues]) |
| 18-Sep-2019 23:19:13 | + conditions['humidity'] = np.mean(relativeHumidity[selectedValues]) |
| 18-Sep-2019 23:19:13 | + if dewPoint is not None: |
| 18-Sep-2019 23:19:13 | + conditions['dewpoint'] = np.mean(dewPoint[selectedValues]) |
| 18-Sep-2019 23:19:13 | + conditions['windspeed'] = np.mean(windSpeed[selectedValues]) |
| 18-Sep-2019 23:19:13 | + conditions['winddirection'] = (180./math.pi)*np.arctan2(np.mean(sinWindDirection[selectedValues]),np.mean(cosWindDirection[selectedValues])) |
| 18-Sep-2019 23:19:13 | + if (conditions['winddirection'] < 0): |
| 18-Sep-2019 23:19:13 | + conditions['winddirection'] += 360 |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Mean weather values for scan %s (field %s)" % (listscan,listfield)) |
| 18-Sep-2019 23:19:13 | + print(" Pressure = %.2f mb" % (conditions['pressure'])) |
| 18-Sep-2019 23:19:13 | + print(" Temperature = %.2f C" % (conditions['temperature'])) |
| 18-Sep-2019 23:19:13 | + if dewPoint is not None: |
| 18-Sep-2019 23:19:13 | + print(" Dew point = %.2f C" % (conditions['dewpoint'])) |
| 18-Sep-2019 23:19:13 | + print(" Relative Humidity = %.2f %%" % (conditions['humidity'])) |
| 18-Sep-2019 23:19:13 | + print(" Wind speed = %.2f m/s" % (conditions['windspeed'])) |
| 18-Sep-2019 23:19:13 | + print(" Wind direction = %.2f deg" % (conditions['winddirection'])) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + return([conditions,myTimes]) |
| 18-Sep-2019 23:19:13 | + # end of getWeather |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def getBasebandNumbers(inputMs) : |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Returns the baseband numbers associated with each spw in the specified ms. |
| 18-Sep-2019 23:19:13 | + Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + if (os.path.exists(inputMs) == False): |
| 18-Sep-2019 23:19:13 | + print("measurement set not found") |
| 18-Sep-2019 23:19:13 | + return -1 |
| 18-Sep-2019 23:19:13 | + mytb = createCasaTool(tbtool) |
| 18-Sep-2019 23:19:13 | + mytb.open("%s/SPECTRAL_WINDOW" % inputMs) |
| 18-Sep-2019 23:19:13 | + if ("BBC_NO" in mytb.colnames()): |
| 18-Sep-2019 23:19:13 | + bbNums = mytb.getcol("BBC_NO") |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + return(-1) |
| 18-Sep-2019 23:19:13 | + mytb.close() |
| 18-Sep-2019 23:19:13 | + return bbNums |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def yigHarmonic(bandString): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Returns the YIG harmonic for the specified ALMA band, given as a string |
| 18-Sep-2019 23:19:13 | + used in casa tables. |
| 18-Sep-2019 23:19:13 | + For example: yigHarmonic('ALMA_RB_03') returns the integer 6. |
| 18-Sep-2019 23:19:13 | + Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + # remove any leading spaces |
| 18-Sep-2019 23:19:13 | + #bandString = bandString[bandString.find('ALMA_RB'):] |
| 18-Sep-2019 23:19:13 | + harmonics = {'ALMA_RB_03':6, 'ALMA_RB_04':6, 'ALMA_RB_06': 18, |
| 18-Sep-2019 23:19:13 | + 'ALMA_RB_07': 18, 'ALMA_RB_08':18, 'ALMA_RB_09':27} |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + harmonic = harmonics[bandString] |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + harmonic = -1 |
| 18-Sep-2019 23:19:13 | + return(harmonic) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def interpretLOs(vis, parentms='', showWVR=False, |
| 18-Sep-2019 23:19:13 | + showCentralFreq=False, verbose=False, show=False, |
| 18-Sep-2019 23:19:13 | + alsoReturnLO2=False, showChannelAverageSpws=False, |
| 18-Sep-2019 23:19:13 | + showOnlyScienceSpws=False, birdieFreq=None, birdieSpw=None, |
| 18-Sep-2019 23:19:13 | + intent='OBSERVE_TARGET#ON_SOURCE', spwsForIntent=None, |
| 18-Sep-2019 23:19:13 | + showEffective=False, showWindowFactors=False, mymsmd=None): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Copied from analysisUtils on May 16, 2017, to replace old version, in order |
| 18-Sep-2019 23:19:13 | + to fix SCOPS-4877. |
| 18-Sep-2019 23:19:13 | + Interpret (and optionally print) the LO settings for an MS from the |
| 18-Sep-2019 23:19:13 | + ASDM_RECEIVER table. |
| 18-Sep-2019 23:19:13 | + Options: |
| 18-Sep-2019 23:19:13 | + showCentralFreq: if True, then show the mean frequency of each spw, |
| 18-Sep-2019 23:19:13 | + otherwise show the frequency of the first channel |
| 18-Sep-2019 23:19:13 | + showWVR: include the WVR spw in the list |
| 18-Sep-2019 23:19:13 | + parentms: if the dataset has been split from a parent dataset, then |
| 18-Sep-2019 23:19:13 | + you may also need to specify the name of the parent ms. |
| 18-Sep-2019 23:19:13 | + alsoReturnLO2: if True, return a second dictionary of the LO2 values |
| 18-Sep-2019 23:19:13 | + birdieFreq: if specified, compute the IF of this RF feature |
| 18-Sep-2019 23:19:13 | + birdieSpw: only necessary if more than one LO1 in the science spws |
| 18-Sep-2019 23:19:13 | + intent: which intent to use in spwsforintent (to find science spws) |
| 18-Sep-2019 23:19:13 | + spwsForIntent: if specified, then avoid the call to spwsforintent |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + Returns: a dictionary of the LO1 values (in Hz) for each spw, keyed by |
| 18-Sep-2019 23:19:13 | + integer. |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + A typical band 7 TDM dataset (prior to splitting) looks like this: |
| 18-Sep-2019 23:19:13 | + SPECTRAL_WINDOW table has 39 rows: row |
| 18-Sep-2019 23:19:13 | + WVR 0 |
| 18-Sep-2019 23:19:13 | + 8 band 3 windows (pointing) 1-8 |
| 18-Sep-2019 23:19:13 | + 8 band 7 windows 9-16 |
| 18-Sep-2019 23:19:13 | + 22 WVR windows 17-38 |
| 18-Sep-2019 23:19:13 | + The corresponding ASDM_RECEIVER table has only 18 rows: |
| 18-Sep-2019 23:19:13 | + WVR 0 |
| 18-Sep-2019 23:19:13 | + 8 band 3 windows 1-8 |
| 18-Sep-2019 23:19:13 | + WVR 9 |
| 18-Sep-2019 23:19:13 | + 8 band 7 windows 10-17 |
| 18-Sep-2019 23:19:13 | + After splitting, the ASDM_RECEIVER table remains the same, but the |
| 18-Sep-2019 23:19:13 | + SPECTRAL WINDOW table then has only 4 rows, as the pointing spws and |
| 18-Sep-2019 23:19:13 | + the channel-averaged data are dropped: |
| 18-Sep-2019 23:19:13 | + 4 band 7 windows |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + lo1s = {} # initialize dictionary to be returned |
| 18-Sep-2019 23:19:13 | + lo2s = {} |
| 18-Sep-2019 23:19:13 | + try: |
| 18-Sep-2019 23:19:13 | + retval = getLOs(vis) |
| 18-Sep-2019 23:19:13 | + [LOs,bands,spws,names,sidebands,receiverIds,spwNames] = retval |
| 18-Sep-2019 23:19:13 | + except: |
| 18-Sep-2019 23:19:13 | + print("getLOs failed") |
| 18-Sep-2019 23:19:13 | + return(retval) |
| 18-Sep-2019 23:19:13 | + if (verbose): print("len(spws) = %d: %s" % (len(spws), str(spws))) |
| 18-Sep-2019 23:19:13 | + maxSpw = np.max(spws) |
| 18-Sep-2019 23:19:13 | + sawWVR = False |
| 18-Sep-2019 23:19:13 | + indices = [] # will exclude the extraneous WVR spws |
| 18-Sep-2019 23:19:13 | + for i in range(len(spws)): |
| 18-Sep-2019 23:19:13 | + if (names[i].find('WVR') >= 0): |
| 18-Sep-2019 23:19:13 | + if (not sawWVR): |
| 18-Sep-2019 23:19:13 | + indices.append(i) |
| 18-Sep-2019 23:19:13 | + sawWVR = True |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + indices.append(i) |
| 18-Sep-2019 23:19:13 | + LOs = np.array(LOs)[indices] |
| 18-Sep-2019 23:19:13 | + bands = np.array(bands)[indices] |
| 18-Sep-2019 23:19:13 | + spws = list(np.array(spws)[indices]) |
| 18-Sep-2019 23:19:13 | + names = np.array(names)[indices] |
| 18-Sep-2019 23:19:13 | + sidebands = np.array(sidebands)[indices] |
| 18-Sep-2019 23:19:13 | + receiverIds = np.array(receiverIds)[indices] |
| 18-Sep-2019 23:19:13 | + index = range(len(spws)) |
| 18-Sep-2019 23:19:13 | + mytb = createCasaTool(tbtool) |
| 18-Sep-2019 23:19:13 | + mytb.open(vis+'/SPECTRAL_WINDOW') |
| 18-Sep-2019 23:19:13 | + # If the data have been split into an ms with fewer spws, then this |
| 18-Sep-2019 23:19:13 | + # table will be smaller (in rows) than the parent MS's table. |
| 18-Sep-2019 23:19:13 | + spwNames = mytb.getcol('NAME') |
| 18-Sep-2019 23:19:13 | + mytb.close() |
| 18-Sep-2019 23:19:13 | + splitted = False |
| 18-Sep-2019 23:19:13 | + if (maxSpw != len(spwNames)-1): |
| 18-Sep-2019 23:19:13 | + splitted = True |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("maxSpw=%d != len(spwNames)=%d)" % (maxSpw, len(spwNames))) |
| 18-Sep-2019 23:19:13 | + if (parentms == '' or parentms == None): |
| 18-Sep-2019 23:19:13 | + print("You appear to have split these data. Please provide the parentms as an argument.") |
| 18-Sep-2019 23:19:13 | + return |
| 18-Sep-2019 23:19:13 | + mytb.open(parentms+'/SPECTRAL_WINDOW') |
| 18-Sep-2019 23:19:13 | + parentSpwNames = mytb.getcol('NAME') |
| 18-Sep-2019 23:19:13 | + mytb.close() |
| 18-Sep-2019 23:19:13 | + extractedRows = [] |
| 18-Sep-2019 23:19:13 | + index = [] |
| 18-Sep-2019 23:19:13 | + for s in range(len(spwNames)): |
| 18-Sep-2019 23:19:13 | + if (len(spwNames[s]) == 0): |
| 18-Sep-2019 23:19:13 | + print("This is an old dataset lacking values in the NAME column of the SPECTRAL_WINDOW table.") |
| 18-Sep-2019 23:19:13 | + return |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("Checking for %s in " % (spwNames[s]), parentSpwNames) |
| 18-Sep-2019 23:19:13 | + extractedRows.append(np.where(parentSpwNames == spwNames[s])[0][0]) |
| 18-Sep-2019 23:19:13 | + index.append(spws.index(extractedRows[-1])) |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("spw %d came from spw %d" % (s, extractedRows[-1])) |
| 18-Sep-2019 23:19:13 | +# extractedRows = the row of the parent SPECTRAL_WINDOW table that matches |
| 18-Sep-2019 23:19:13 | +# the split-out spw |
| 18-Sep-2019 23:19:13 | +# index = the row of the ASDM_RECEIVER table that matches the split-out spw |
| 18-Sep-2019 23:19:13 | + vis = parentms |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("spwNames = ", spwNames) |
| 18-Sep-2019 23:19:13 | + print("spws = ", spws) |
| 18-Sep-2019 23:19:13 | + print("bands = ", bands) |
| 18-Sep-2019 23:19:13 | + output = "LOs = " |
| 18-Sep-2019 23:19:13 | + for LO in LOs: |
| 18-Sep-2019 23:19:13 | + output += "%.3f, " % (LO[0]*1e-9) |
| 18-Sep-2019 23:19:13 | + print(output) |
| 18-Sep-2019 23:19:13 | + print("names = ", names) |
| 18-Sep-2019 23:19:13 | + print("index = ", index) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + bbc = getBasebandNumbers(vis) # does not use msmd |
| 18-Sep-2019 23:19:13 | + if (show): |
| 18-Sep-2019 23:19:13 | + print('Row refers to the row number in the ASDM_RECEIVER table (starting at 0).') |
| 18-Sep-2019 23:19:13 | + if (showCentralFreq): |
| 18-Sep-2019 23:19:13 | + myline = 'Row spw BB RxBand CenFreq Nchan LO1(GHz) LO2(GHz) Sampler YIG(GHz) TFBoffset(MHz)' |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + myline = 'Row spw BB RxBand Ch1Freq Nchan LO1(GHz) LO2(GHz) Sampler YIG(GHz) TFBoffset(MHz)' |
| 18-Sep-2019 23:19:13 | + if (showEffective): |
| 18-Sep-2019 23:19:13 | + myline += ' Eff.BW(MHz) Res(MHz) Width(MHz)' |
| 18-Sep-2019 23:19:13 | + if (showWindowFactors): |
| 18-Sep-2019 23:19:13 | + myline += ' windowFactors' |
| 18-Sep-2019 23:19:13 | + print(myline) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + # Loop over all rows in the ASDM_RECEIVER table, unless we've split, in |
| 18-Sep-2019 23:19:13 | + # which case this will loop over the N spws in the table. |
| 18-Sep-2019 23:19:13 | + needToClose = False |
| 18-Sep-2019 23:19:13 | + if mymsmd is None or mymsmd == '': |
| 18-Sep-2019 23:19:13 | + mymsmd = createCasaTool(msmdtool) |
| 18-Sep-2019 23:19:13 | + mymsmd.open(vis) |
| 18-Sep-2019 23:19:13 | + needToClose = True |
| 18-Sep-2019 23:19:13 | + if (spwsForIntent == None): |
| 18-Sep-2019 23:19:13 | + if intent in mymsmd.intents(): # prevent a warning of OBSERVE_TARGET does not exist |
| 18-Sep-2019 23:19:13 | + scienceSpws = np.setdiff1d(mymsmd.spwsforintent(intent),mymsmd.wvrspws()) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + scienceSpws = [] |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + scienceSpws = spwsForIntent |
| 18-Sep-2019 23:19:13 | + birdieIF = 0 |
| 18-Sep-2019 23:19:13 | + if (birdieFreq is not None): |
| 18-Sep-2019 23:19:13 | + birdieFreq = parseFrequencyArgumentToHz(birdieFreq) |
| 18-Sep-2019 23:19:13 | + birdieFreqGHz = parseFrequencyArgumentToGHz(birdieFreq) |
| 18-Sep-2019 23:19:13 | + fdmSpws = mymsmd.almaspws(fdm=True) |
| 18-Sep-2019 23:19:13 | + for i in range(len(index)): |
| 18-Sep-2019 23:19:13 | + if (verbose): |
| 18-Sep-2019 23:19:13 | + print("index[%d]=%d" % (i,index[i])) |
| 18-Sep-2019 23:19:13 | + print("spws[%d] = %d" % (index[i], spws[index[i]])) |
| 18-Sep-2019 23:19:13 | + myspw = i |
| 18-Sep-2019 23:19:13 | + if (birdieFreq is not None and birdieIF == 0): |
| 18-Sep-2019 23:19:13 | + if (myspw == birdieSpw): |
| 18-Sep-2019 23:19:13 | + if verbose: |
| 18-Sep-2019 23:19:13 | + print("spw=%d, Computing IF = %f - %f" % (myspw, birdieFreq, LOs[index[i]][0])) |
| 18-Sep-2019 23:19:13 | + birdieIF = np.fabs(birdieFreq - LOs[index[i]][0]) |
| 18-Sep-2019 23:19:13 | + elif (myspw in scienceSpws and birdieSpw==None): |
| 18-Sep-2019 23:19:13 | + if verbose: |
| 18-Sep-2019 23:19:13 | + print("spw=%d (in %s), Computing IF = %f - %f" % (myspw, str(scienceSpws), birdieFreq, LOs[index[i]][0])) |
| 18-Sep-2019 23:19:13 | + birdieIF = np.fabs(birdieFreq - LOs[index[i]][0]) |
| 18-Sep-2019 23:19:13 | + freqs = mymsmd.chanfreqs(myspw) |
| 18-Sep-2019 23:19:13 | + meanFreqGHz = mymsmd.meanfreq(myspw) * (1e-9) |
| 18-Sep-2019 23:19:13 | + if (myspw not in scienceSpws and showOnlyScienceSpws): continue |
| 18-Sep-2019 23:19:13 | + if (len(freqs) < 2 and showChannelAverageSpws==False): |
| 18-Sep-2019 23:19:13 | + continue |
| 18-Sep-2019 23:19:13 | + if (bands[index[i]].split('_')[-1].isdigit()): |
| 18-Sep-2019 23:19:13 | + rxband = bands[index[i]].split('_')[-1] |
| 18-Sep-2019 23:19:13 | + elif (showWVR): |
| 18-Sep-2019 23:19:13 | + rxband = 'WVR' |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + continue |
| 18-Sep-2019 23:19:13 | + line = "%2d %2d %d %3s " % (spws[index[i]], myspw, bbc[myspw], rxband) |
| 18-Sep-2019 23:19:13 | + if (showCentralFreq): |
| 18-Sep-2019 23:19:13 | + line += "%10.6f %4d " % (meanFreqGHz,len(freqs)) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + line += "%10.6f %4d " % (freqs[0],len(freqs)) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + if (LOs[index[i]][0] < 0): |
| 18-Sep-2019 23:19:13 | + print(line) |
| 18-Sep-2019 23:19:13 | + continue |
| 18-Sep-2019 23:19:13 | + if (bbc[myspw] > 0): |
| 18-Sep-2019 23:19:13 | + if (splitted): |
| 18-Sep-2019 23:19:13 | + lo1s[i] = LOs[index[i]][0] |
| 18-Sep-2019 23:19:13 | + lo2s[i] = LOs[index[i]][1] |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + lo1s[myspw] = LOs[index[i]][0] |
| 18-Sep-2019 23:19:13 | + lo2s[myspw] = LOs[index[i]][1] |
| 18-Sep-2019 23:19:13 | + for j in range(len(LOs[index[i]])): |
| 18-Sep-2019 23:19:13 | + if (j != 2): |
| 18-Sep-2019 23:19:13 | + line = line + '%10.6f' % (LOs[index[i]][j]*1e-9) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + line = line + '%5.2f' % (LOs[index[i]][j]*1e-9) |
| 18-Sep-2019 23:19:13 | + yig = LOs[index[i]][0] / yigHarmonic(bands[index[i]]) |
| 18-Sep-2019 23:19:13 | + if (yig > 0): |
| 18-Sep-2019 23:19:13 | + line = line + ' %.6f' % (yig*1e-9) |
| 18-Sep-2019 23:19:13 | + if (myspw in fdmSpws): |
| 18-Sep-2019 23:19:13 | + # work out what LO4 must have been |
| 18-Sep-2019 23:19:13 | + LO1 = LOs[index[i]][0] |
| 18-Sep-2019 23:19:13 | + LO2 = LOs[index[i]][1] |
| 18-Sep-2019 23:19:13 | + LO3 = LOs[index[i]][2] |
| 18-Sep-2019 23:19:13 | + if (sidebands[index[i]][0] == 'USB'): |
| 18-Sep-2019 23:19:13 | + IFlocation = LO3 - (LO2 - (meanFreqGHz*1e9 - LO1)) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + IFlocation = LO3 - (LO2 - (LO1 - meanFreqGHz*1e9)) |
| 18-Sep-2019 23:19:13 | + LO4 = 2e9 + IFlocation |
| 18-Sep-2019 23:19:13 | + TFBLOoffset = LO4-3e9 |
| 18-Sep-2019 23:19:13 | + line += '%9.3f %+8.3f ' % (LO4 * 1e-6, TFBLOoffset * 1e-6) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + line += 19*' ' |
| 18-Sep-2019 23:19:13 | + if (showEffective): |
| 18-Sep-2019 23:19:13 | + line += '%9.4f %9.4f %9.4f' % (effectiveBandwidth(vis, myspw)*1e-6, |
| 18-Sep-2019 23:19:13 | + effectiveResolution(vis, myspw)*1e-6, |
| 18-Sep-2019 23:19:13 | + getChanWidths(mymsmd, myspw)*1e-6) |
| 18-Sep-2019 23:19:13 | + if (showWindowFactors): |
| 18-Sep-2019 23:19:13 | + chanwidth = abs(getChanWidths(mymsmd,myspw)) |
| 18-Sep-2019 23:19:13 | + line += ' %.4f %.4f' % (effectiveBandwidth(vis, myspw)/chanwidth, |
| 18-Sep-2019 23:19:13 | + effectiveResolution(vis, myspw)/chanwidth) |
| 18-Sep-2019 23:19:13 | + if (bands[index[i]].find('ALMA_RB_06')>=0 or bands[index[i]].find('ALMA_RB_09')>=0): |
| 18-Sep-2019 23:19:13 | + if (len(LOs[index[i]]) > 1): |
| 18-Sep-2019 23:19:13 | + if (LOs[index[i]][1] < 11.3e9 and LOs[index[i]][1] > 10.5e9): |
| 18-Sep-2019 23:19:13 | + line = line + ' leakage of LO2 undesired sideband may degrade dynamic range' |
| 18-Sep-2019 23:19:13 | + if (bands[index[i]].find('ALMA_RB_06')>=0): |
| 18-Sep-2019 23:19:13 | + line += ' (and YIG may leak in)' |
| 18-Sep-2019 23:19:13 | + yigLeakage = LOs[index[i]][0] + (LOs[index[i]][1] - LOs[index[i]][2]) + (yig - LOs[index[i]][1]) |
| 18-Sep-2019 23:19:13 | + if (yigLeakage > 0): |
| 18-Sep-2019 23:19:13 | + line = line + ' at %.6f' % (yigLeakage*1e-9) |
| 18-Sep-2019 23:19:13 | + if (show): print(line) |
| 18-Sep-2019 23:19:13 | + if needToClose: |
| 18-Sep-2019 23:19:13 | + mymsmd.done() |
| 18-Sep-2019 23:19:13 | + if (birdieIF != 0): |
| 18-Sep-2019 23:19:13 | + print("The feature at %f GHz is at IF = %f GHz." % (birdieFreqGHz, birdieIF*1e-9)) |
| 18-Sep-2019 23:19:13 | + if (alsoReturnLO2): |
| 18-Sep-2019 23:19:13 | + return(lo1s, lo2s) |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + return(lo1s) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def mjdSecondsToMJDandUT(mjdsec): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Converts a value of MJD seconds into MJD, and into a UT date/time string. |
| 18-Sep-2019 23:19:13 | + example: (56000.0, '2012-03-14 00:00:00 UT') |
| 18-Sep-2019 23:19:13 | + Caveat: only works for a scalar input value |
| 18-Sep-2019 23:19:13 | + Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + myme = createCasaTool(metool) |
| 18-Sep-2019 23:19:13 | + today = myme.epoch('utc','today') |
| 18-Sep-2019 23:19:13 | + mjd = np.array(mjdsec) / 86400. |
| 18-Sep-2019 23:19:13 | + today['m0']['value'] = mjd |
| 18-Sep-2019 23:19:13 | + myqa = createCasaTool(qatool) |
| 18-Sep-2019 23:19:13 | + hhmmss = myqa.time(today['m0'], form='', prec=0, showform=False)[0] |
| 18-Sep-2019 23:19:13 | +# print("hhmmss = ", hhmmss) |
| 18-Sep-2019 23:19:13 | + date = myqa.splitdate(today['m0']) |
| 18-Sep-2019 23:19:13 | + myqa.done() |
| 18-Sep-2019 23:19:13 | + utstring = "%s-%02d-%02d %s UT" % (date['year'],date['month'],date['monthday'],hhmmss) |
| 18-Sep-2019 23:19:13 | + myme.done() |
| 18-Sep-2019 23:19:13 | + return(mjd, utstring) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def ComputeSolarAzElForObservatory(mjdsec, mymsmd): |
| 18-Sep-2019 23:19:13 | + pos = mymsmd.observatoryposition() |
| 18-Sep-2019 23:19:13 | + longitude = pos['m0']['value'] * 180/np.pi |
| 18-Sep-2019 23:19:13 | + latitude = pos['m1']['value'] * 180/np.pi |
| 18-Sep-2019 23:19:13 | + return(ComputeSolarAzElLatLong(mjdsec,latitude,longitude)) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def ComputeSolarAzElLatLong(mjdsec,latitude,longitude): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Computes the apparent Az,El of the Sun for a specified time and location |
| 18-Sep-2019 23:19:13 | + on Earth. Latitude and longitude must arrive in degrees, with positive |
| 18-Sep-2019 23:19:13 | + longitude meaning east of Greenwich. |
| 18-Sep-2019 23:19:13 | + -- Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + DEG_TO_RAD = math.pi/180. |
| 18-Sep-2019 23:19:13 | + RAD_TO_DEG = 180/math.pi |
| 18-Sep-2019 23:19:13 | + HRS_TO_RAD = math.pi/12. |
| 18-Sep-2019 23:19:13 | + [RA,Dec] = ComputeSolarRADec(mjdsec) |
| 18-Sep-2019 23:19:13 | + LST = ComputeLST(mjdsec, longitude) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + phi = latitude*DEG_TO_RAD |
| 18-Sep-2019 23:19:13 | + hourAngle = HRS_TO_RAD*(LST - RA) |
| 18-Sep-2019 23:19:13 | + azimuth = RAD_TO_DEG*math.atan2(math.sin(hourAngle), (math.cos(hourAngle)*math.sin(phi) - math.tan(Dec*DEG_TO_RAD)*math.cos(phi))) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + # the following is to convert from South=0 (which the French formula uses) |
| 18-Sep-2019 23:19:13 | + # to North=0, which is what the rest of the world uses */ |
| 18-Sep-2019 23:19:13 | + azimuth += 180.0; |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + if (azimuth > 360.0): |
| 18-Sep-2019 23:19:13 | + azimuth -= 360.0 |
| 18-Sep-2019 23:19:13 | + if (azimuth < 0.0): |
| 18-Sep-2019 23:19:13 | + azimuth += 360.0 |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + argument = math.sin(phi)*math.sin(Dec*DEG_TO_RAD) + math.cos(phi)*math.cos(Dec*DEG_TO_RAD) * math.cos(hourAngle); |
| 18-Sep-2019 23:19:13 | + elevation = RAD_TO_DEG*math.asin(argument); |
| 18-Sep-2019 23:19:13 | + return([azimuth,elevation]) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def ComputeSolarRADec(mjdsec): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Computes the RA,Dec of the Sun for a specified time. -- Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + jd = mjdToJD(mjdsec/86400.) |
| 18-Sep-2019 23:19:13 | + RAD_TO_DEG = 180/math.pi |
| 18-Sep-2019 23:19:13 | + RAD_TO_HRS = (1.0/0.2617993877991509) |
| 18-Sep-2019 23:19:13 | + DEG_TO_RAD = math.pi/180. |
| 18-Sep-2019 23:19:13 | + T = (jd - 2451545.0) / 36525.0 |
| 18-Sep-2019 23:19:13 | + Lo = 280.46646 + 36000.76983*T + 0.0003032*T*T |
| 18-Sep-2019 23:19:13 | + M = 357.52911 + 35999.05029*T - 0.0001537*T*T |
| 18-Sep-2019 23:19:13 | + Mrad = M * DEG_TO_RAD |
| 18-Sep-2019 23:19:13 | + e = 0.016708634 - 0.000042037*T - 0.0000001267*T*T |
| 18-Sep-2019 23:19:13 | + C = (1.914602 - 0.004817*T - 0.000014*T*T) * math.sin(Mrad) + (0.019993 - 0.000101*T) * math.sin(2*Mrad) + 0.000289*math.sin(3*Mrad) |
| 18-Sep-2019 23:19:13 | + L = Lo + C |
| 18-Sep-2019 23:19:13 | + nu = DEG_TO_RAD*(M + C) |
| 18-Sep-2019 23:19:13 | + R = 1.000001018 * (1-e*e) / (1 + e*math.cos(nu)) |
| 18-Sep-2019 23:19:13 | + Omega = DEG_TO_RAD*(125.04 - 1934.136*T) |
| 18-Sep-2019 23:19:13 | + mylambda = DEG_TO_RAD*(L - 0.00569 - 0.00478 * math.sin(Omega)) |
| 18-Sep-2019 23:19:13 | + epsilon0 = (84381.448 - 46.8150*T - 0.00059*T*T + 0.001813*T*T*T) / 3600. |
| 18-Sep-2019 23:19:13 | + epsilon = (epsilon0 + 0.00256 * math.cos(Omega)) * DEG_TO_RAD |
| 18-Sep-2019 23:19:13 | + rightAscension = RAD_TO_HRS*math.atan2(math.cos(epsilon)*math.sin(mylambda), math.cos(mylambda)) |
| 18-Sep-2019 23:19:13 | + if (rightAscension < 0): |
| 18-Sep-2019 23:19:13 | + rightAscension += 24.0 |
| 18-Sep-2019 23:19:13 | + argument = math.sin(epsilon) * math.sin(mylambda) |
| 18-Sep-2019 23:19:13 | + declination = RAD_TO_DEG*math.asin(argument) |
| 18-Sep-2019 23:19:13 | + return([rightAscension, declination]) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def angularSeparation(ra0,dec0,ra1,dec1, returnComponents=False): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Usage: angularSeparation(ra0,dec0,ra1,dec1) |
| 18-Sep-2019 23:19:13 | + Computes the great circle angle between two celestial coordinates. |
| 18-Sep-2019 23:19:13 | + using the Vincenty formula (from wikipedia) which is correct for all |
| 18-Sep-2019 23:19:13 | + angles, as long as you use atan2() to handle a zero denominator. |
| 18-Sep-2019 23:19:13 | + See http://en.wikipedia.org/wiki/Great_circle_distance |
| 18-Sep-2019 23:19:13 | + ra,dec must be given in degrees, as is the output. |
| 18-Sep-2019 23:19:13 | + It also works for the az,el coordinate system. |
| 18-Sep-2019 23:19:13 | + Comopnent separations are field_0 minus field_1. |
| 18-Sep-2019 23:19:13 | + See also angularSeparationRadians() |
| 18-Sep-2019 23:19:13 | + -- Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + ra0 *= math.pi/180. |
| 18-Sep-2019 23:19:13 | + dec0 *= math.pi/180. |
| 18-Sep-2019 23:19:13 | + ra1 *= math.pi/180. |
| 18-Sep-2019 23:19:13 | + dec1 *= math.pi/180. |
| 18-Sep-2019 23:19:13 | + deltaLong = ra0-ra1 |
| 18-Sep-2019 23:19:13 | + argument1 = (((math.cos(dec1)*math.sin(deltaLong))**2) + |
| 18-Sep-2019 23:19:13 | + ((math.cos(dec0)*math.sin(dec1)-math.sin(dec0)*math.cos(dec1)*math.cos(deltaLong))**2))**0.5 |
| 18-Sep-2019 23:19:13 | + argument2 = math.sin(dec0)*math.sin(dec1) + math.cos(dec0)*math.cos(dec1)*math.cos(deltaLong) |
| 18-Sep-2019 23:19:13 | + angle = math.atan2(argument1, argument2) / (math.pi/180.) |
| 18-Sep-2019 23:19:13 | + if (returnComponents): |
| 18-Sep-2019 23:19:13 | + radegrees = (ra0-ra1)*180/math.pi |
| 18-Sep-2019 23:19:13 | + decdegrees = (dec0-dec1)*180/math.pi |
| 18-Sep-2019 23:19:13 | + retval = angle,radegrees,decdegrees |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + retval = angle |
| 18-Sep-2019 23:19:13 | + return(retval) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def ComputeDewPointCFromRHAndTempC(relativeHumidity, temperature): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + inputs: relativeHumidity in percentage, temperature in C |
| 18-Sep-2019 23:19:13 | + output: in degrees C |
| 18-Sep-2019 23:19:13 | + Uses formula from http://en.wikipedia.org/wiki/Dew_point#Calculating_the_dew_point |
| 18-Sep-2019 23:19:13 | + Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + temperature = np.array(temperature) # protect against it being a list |
| 18-Sep-2019 23:19:13 | + relativeHumidity = np.array(relativeHumidity) # protect against it being a list |
| 18-Sep-2019 23:19:13 | + es = 6.112*np.exp(17.67*temperature/(temperature+243.5)) |
| 18-Sep-2019 23:19:13 | + E = relativeHumidity*0.01*es |
| 18-Sep-2019 23:19:13 | + dewPoint = 243.5*np.log(E/6.112)/(17.67-np.log(E/6.112)) |
| 18-Sep-2019 23:19:13 | + return(dewPoint) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def ComputeLST(mjdsec, longitude): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Computes the LST (in hours) for a specified time and longitude. |
| 18-Sep-2019 23:19:13 | + The input longitude is in degrees, where east of Greenwich is positive. |
| 18-Sep-2019 23:19:13 | + -- Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + JD = mjdToJD(mjdsec/86400.) |
| 18-Sep-2019 23:19:13 | + T = (JD - 2451545.0) / 36525.0 |
| 18-Sep-2019 23:19:13 | + sidereal = 280.46061837 + 360.98564736629*(JD - 2451545.0) + 0.000387933*T*T - T*T*T/38710000. |
| 18-Sep-2019 23:19:13 | + # now we have LST in Greenwich, need to scale back to site |
| 18-Sep-2019 23:19:13 | + sidereal += longitude |
| 18-Sep-2019 23:19:13 | + sidereal /= 360. |
| 18-Sep-2019 23:19:13 | + sidereal -= np.floor(sidereal) |
| 18-Sep-2019 23:19:13 | + sidereal *= 24.0 |
| 18-Sep-2019 23:19:13 | + if (sidereal < 0): |
| 18-Sep-2019 23:19:13 | + sidereal += 24 |
| 18-Sep-2019 23:19:13 | + if (sidereal >= 24): |
| 18-Sep-2019 23:19:13 | + sidereal -= 24 |
| 18-Sep-2019 23:19:13 | + return(sidereal) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def mjdToJD(MJD): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Converts an MJD value to JD |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + JD = MJD + 2400000.5 |
| 18-Sep-2019 23:19:13 | + return(JD) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def splitListIntoContiguousLists(mylist): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Converts [1,2,3,5,6,7] into [[1,2,3],[5,6,7]], etc. |
| 18-Sep-2019 23:19:13 | + -Todd Hunter |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + mylists = [] |
| 18-Sep-2019 23:19:13 | + newlist = [mylist[0]] |
| 18-Sep-2019 23:19:13 | + for i in range(1,len(mylist)): |
| 18-Sep-2019 23:19:13 | + if (mylist[i-1] != mylist[i]-1): |
| 18-Sep-2019 23:19:13 | + mylists.append(newlist) |
| 18-Sep-2019 23:19:13 | + newlist = [mylist[i]] |
| 18-Sep-2019 23:19:13 | + else: |
| 18-Sep-2019 23:19:13 | + newlist.append(mylist[i]) |
| 18-Sep-2019 23:19:13 | + mylists.append(newlist) |
| 18-Sep-2019 23:19:13 | + return(mylists) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +# Removed for CAS-8065 |
| 18-Sep-2019 23:19:13 | +#def getScansForTimes(mymsmd, scantimes): |
| 18-Sep-2019 23:19:13 | +# myscans = [] |
| 18-Sep-2019 23:19:13 | +# myscantimes = [] |
| 18-Sep-2019 23:19:13 | +## print("len(scantimes) = ", len(scantimes)) |
| 18-Sep-2019 23:19:13 | +# scantimes = splitListIntoContiguousLists(scantimes) |
| 18-Sep-2019 23:19:13 | +# for t in scantimes: |
| 18-Sep-2019 23:19:13 | +# mean_t = np.mean(t) |
| 18-Sep-2019 23:19:13 | +# range_t = (1+np.max(t)-np.min(t))*0.5 |
| 18-Sep-2019 23:19:13 | +# scans_t = mymsmd.scansfortimes(mean_t, range_t) |
| 18-Sep-2019 23:19:13 | +# if (len(scans_t) > 0): |
| 18-Sep-2019 23:19:13 | +# scan = scans_t[0] |
| 18-Sep-2019 23:19:13 | +# # print("scansfortime(%f) = " % (t), scan) |
| 18-Sep-2019 23:19:13 | +# myscans.append(scan) |
| 18-Sep-2019 23:19:13 | +# myscantimes.append(t) |
| 18-Sep-2019 23:19:13 | +# return(myscans, myscantimes) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def pruneFilelist(filelist): |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + Reduce size of filenames in filelist to the extent that current working directory |
| 18-Sep-2019 23:19:13 | + agrees with the path. |
| 18-Sep-2019 23:19:13 | + """ |
| 18-Sep-2019 23:19:13 | + mypwd = os.getcwd() + '/' |
| 18-Sep-2019 23:19:13 | + newfilelist = [] |
| 18-Sep-2019 23:19:13 | + for f in filelist: |
| 18-Sep-2019 23:19:13 | + fstart = 0 |
| 18-Sep-2019 23:19:13 | + if (f.find(mypwd) == 0): |
| 18-Sep-2019 23:19:13 | + fstart = len(mypwd) |
| 18-Sep-2019 23:19:13 | + newfilelist.append(f[fstart:]) |
| 18-Sep-2019 23:19:13 | + return(newfilelist) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | diff --git a/tests/tasks/test_fringefit.py b/tests/tasks/test_fringefit.py |
| 18-Sep-2019 23:19:13 | new file mode 100644 |
| 18-Sep-2019 23:19:13 | index 0000000..e0c4186 |
| 18-Sep-2019 23:19:13 | --- /dev/null |
| 18-Sep-2019 23:19:13 | +++ b/tests/tasks/test_fringefit.py |
| 18-Sep-2019 23:19:13 | @@ -0,0 +1,67 @@ |
| 18-Sep-2019 23:19:13 | +from __future__ import absolute_import |
| 18-Sep-2019 23:19:13 | +from __future__ import print_function |
| 18-Sep-2019 23:19:13 | +import os |
| 18-Sep-2019 23:19:13 | +import sys |
| 18-Sep-2019 23:19:13 | +import shutil |
| 18-Sep-2019 23:19:13 | +import unittest |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +# is_CASA6 and is_python3 |
| 18-Sep-2019 23:19:13 | +from casatasks.private.casa_transition import * |
| 18-Sep-2019 23:19:13 | +if is_CASA6: |
| 18-Sep-2019 23:19:13 | + from casatools import ms, ctsys |
| 18-Sep-2019 23:19:13 | + from casatasks import fringefit |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + ctsys_resolve = ctsys.resolve |
| 18-Sep-2019 23:19:13 | +else: |
| 18-Sep-2019 23:19:13 | + from __main__ import default |
| 18-Sep-2019 23:19:13 | + from tasks import * |
| 18-Sep-2019 23:19:13 | + from taskinit import * |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + dataRoot = os.path.join(os.environ.get('CASAPATH').split()[0],'data') |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + def ctsys_resolve(apath): |
| 18-Sep-2019 23:19:13 | + return os.path.join(dataRoot,apath) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +if is_python3: |
| 18-Sep-2019 23:19:13 | + ### for testhelper import |
| 18-Sep-2019 23:19:13 | + sys.path.append(os.path.abspath(os.path.dirname(__file__))) |
| 18-Sep-2019 23:19:13 | + import testhelper as th |
| 18-Sep-2019 23:19:13 | +else: |
| 18-Sep-2019 23:19:13 | + import testhelper as th |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +datapath = ctsys_resolve('regression/evn') |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +class Fringefit_tests(unittest.TestCase): |
| 18-Sep-2019 23:19:13 | + prefix = 'n08c1' |
| 18-Sep-2019 23:19:13 | + msfile = prefix + '.ms' |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + def setUp(self): |
| 18-Sep-2019 23:19:13 | + shutil.copytree(os.path.join(datapath, self.msfile), self.msfile) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + def tearDown(self): |
| 18-Sep-2019 23:19:13 | + shutil.rmtree(self.msfile) |
| 18-Sep-2019 23:19:13 | + shutil.rmtree(self.prefix + '.sbdcal', True) |
| 18-Sep-2019 23:19:13 | + shutil.rmtree(self.prefix + '-zerorates.sbdcal', True) |
| 18-Sep-2019 23:19:13 | + shutil.rmtree(self.prefix + '.mbdcal', True) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + def test_sbd(self): |
| 18-Sep-2019 23:19:13 | + sbdcal = self.prefix + '.sbdcal' |
| 18-Sep-2019 23:19:13 | + fringefit(vis=self.msfile, caltable=sbdcal, refant='EF') |
| 18-Sep-2019 23:19:13 | + reference = os.path.join(datapath, sbdcal) |
| 18-Sep-2019 23:19:13 | + self.assertTrue(th.compTables(sbdcal, reference, ['WEIGHT'])) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | + def test_mbd(self): |
| 18-Sep-2019 23:19:13 | + sbdcal = self.prefix + '-zerorates.sbdcal' |
| 18-Sep-2019 23:19:13 | + mbdcal = self.prefix + '.mbdcal' |
| 18-Sep-2019 23:19:13 | + fringefit(vis=self.msfile, caltable=sbdcal, field='4C39.25', |
| 18-Sep-2019 23:19:13 | + refant='EF', zerorates=True) |
| 18-Sep-2019 23:19:13 | + fringefit(vis=self.msfile, caltable=mbdcal, field='J0916+3854', |
| 18-Sep-2019 23:19:13 | + combine='spw', gaintable=[sbdcal], refant='EF') |
| 18-Sep-2019 23:19:13 | + reference = os.path.join(datapath, mbdcal) |
| 18-Sep-2019 23:19:13 | + self.assertTrue(th.compTables(mbdcal, reference, ['WEIGHT'])) |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +def suite(): |
| 18-Sep-2019 23:19:13 | + return [Fringefit_tests] |
| 18-Sep-2019 23:19:13 | + |
| 18-Sep-2019 23:19:13 | +if __name__ == '__main__': |
| 18-Sep-2019 23:19:13 | + unittest.main() |
| 18-Sep-2019 23:19:14 | Argument values: |
| 18-Sep-2019 23:19:14 | casatasks |
| 18-Sep-2019 23:19:14 | /export/home/cbt-el7-6/bamboohome/xml-data/build-dir/CASA-CTB-JOB1/casatasks |
| 18-Sep-2019 23:19:14 | master |
| 18-Sep-2019 23:19:14 | None |
| 18-Sep-2019 23:19:14 | Last Tag: 2019.150 |
| 18-Sep-2019 23:19:14 | New Tag: 2019.151 |
| 18-Sep-2019 23:19:14 | Creating a new tag |
| 18-Sep-2019 23:19:14 | |
| 18-Sep-2019 23:19:14 | None |
| 18-Sep-2019 23:19:14 | Push enabled. Pushing new tag upstream. |
| 18-Sep-2019 23:19:14 | To ssh://git@open-bitbucket.nrao.edu:7999/casa/casatasks.git |
| 18-Sep-2019 23:19:14 | * [new tag] 2019.151 -> 2019.151 |
| 18-Sep-2019 23:19:14 | |
| 18-Sep-2019 23:19:14 | None |
| 18-Sep-2019 23:19:14 | Finished task 'Create wheel' with result: Success |
| 18-Sep-2019 23:19:14 | Running post build plugin 'Docker Container Cleanup' |
| 18-Sep-2019 23:19:14 | Beginning to execute external process for build 'CASA - Casa6 Tasks Wheel Build - Build and Package RHEL6 #146 (CASA-CTB-JOB1-146)' ... running command line: /bin/docker rm -f task-wheel-container ... in: /export/home/cbt-el7-6/bamboohome/xml-data/build-dir/CASA-CTB-JOB1 |
| 18-Sep-2019 23:19:15 | task-wheel-container |
| 18-Sep-2019 23:19:15 | Running post build plugin 'NCover Results Collector' |
| 18-Sep-2019 23:19:15 | Running post build plugin 'Clover Results Collector' |
| 18-Sep-2019 23:19:15 | Running post build plugin 'npm Cache Cleanup' |
| 18-Sep-2019 23:19:15 | Running post build plugin 'Artifact Copier' |
| 18-Sep-2019 23:19:15 | Publishing an artifact: EL6 wheel |
| 18-Sep-2019 23:19:15 | Finished publishing of artifact Non required shared artifact: [EL6 wheel], pattern: [*-any.whl] anchored at: [wheeldirectory/] in 165.7 ms |
| 18-Sep-2019 23:19:15 | Finalising the build... |
| 18-Sep-2019 23:19:15 | Stopping timer. |
| 18-Sep-2019 23:19:15 | Build CASA-CTB-JOB1-146 completed. |
| 18-Sep-2019 23:19:15 | Running on server: post build plugin 'NCover Results Collector' |
| 18-Sep-2019 23:19:15 | Running on server: post build plugin 'Build Hanging Detection Configuration' |
| 18-Sep-2019 23:19:15 | Running on server: post build plugin 'Clover Delta Calculator' |
| 18-Sep-2019 23:19:15 | Running on server: post build plugin 'Maven Dependencies Postprocessor' |
| 18-Sep-2019 23:19:15 | All post build plugins have finished |
| 18-Sep-2019 23:19:15 | Generating build results summary... |
| 18-Sep-2019 23:19:15 | Saving build results to disk... |
| 18-Sep-2019 23:19:15 | Logging substituted variables... |
| 18-Sep-2019 23:19:15 | Indexing build results... |
| 18-Sep-2019 23:19:15 | Finished building CASA-CTB-JOB1-146. |