# Albus_Iono.py # Python stuff for dealing with Ionosphere information in AIPS AI tables # 2006 Jan 05 James M Anderson --JIVE start ################################################################################ # some import commands The User should not need to change this. ################################################################################ ################################################################################ # Mark's ParselTongue stuff from AIPS import AIPS, AIPSDisk from AIPSTask import AIPSTask, AIPSList from AIPSData import AIPSUVData, AIPSImage from AIPSTV import AIPSTV import Wizardry.AIPSData ################################################################################ # miscellaneous stuff import copy, optparse, os, sys import re, string import numarray import numarray.ma import numarray.ieeespecial import numarray.nd_image import inspect import warnings import math import time as systime ################################################################################ # JMA's ionosphere stuff import AlbusIonosphere #import AIPSUVData, AIPSTableRow ################################################################################ # Global variables SECONDS_PER_DAY = 24.0*60.0*60.0 DAYS_PER_SECOND = 1.0/SECONDS_PER_DAY M_DEG2RAD = math.pi/180.0 M_RAD2DEG = 180.0/math.pi # See JMA Notes from 2006 Jan 10 and some of BOb Campbell's notes SPEED_OF_LIGHT = 299792458.0 # m s^{-1} IONOSPHERE_Kp_2 = 80.6163848 # m^3 s^{-2} IONOSPHERE_KB = 2.79924925E10 # s^{-1} T^{-1} # AIPS defines things at \lambda = 1.0 m # TEC in units of electrons m^{-2} # RM in units of electrons m^{-2} T AIPS_TEC_TO_DISP = IONOSPHERE_Kp_2 / (2.0 * SPEED_OF_LIGHT**3) AIPS_RM_TO_RM = math.pi * IONOSPHERE_Kp_2 * IONOSPHERE_KB / (SPEED_OF_LIGHT**3) ################################################################################ # Get the day of year from the Year, month, day def get_day_of_year(year, month, day): """Get the day-of-year integer from the year/month/day year I YYYY month I MM starting from January == 1 day I DD starting from 1 == 1 """ day_of_year_list = [0,31,59,90,120,151,181,212,243,273,304,334] doy = day_of_year_list[month-1] + day if(month>2): if((year&0x3)==0): if((year % 100 != 0) or (year % 400 == 0)): doy = doy+1 return doy ################################################################################ # Get the day of year from the Year, month, day for the start of observations def get_observation_year_month_day(aips_data): date_string = aips_data.header.date_obs date_list = date_string.split('-') year = int(date_list[0]) month = int(date_list[1]) day = int(date_list[2]) return (year, month, day) ################################################################################ # Get the day of year from the Year, month, day for the start of observations def get_observation_day_of_year(aips_data): (year, month, day) = get_observation_year_month_day(aips_data) return get_day_of_year(year, month, day) ################################################################################ # Get the number of days of observations def get_number_days_observations(aips_data): """Find out how long the observations are. Assumes the data has been sorted """ nx_table = aips_data.table('AIPS NX', 0) num_rows = len(nx_table) # The NX table is FORTRANish # Check the first observation first_row = nx_table[0] start_time = first_row['TIME'][0] - 0.5 * first_row['TIME INTERVAL'][0] if(start_time < 0.0): raise RuntimeError, "Starting observation before starting date: Run FXTIM!" # Now check the last scan last_row = nx_table[num_rows-1] end_time = last_row['TIME'][0] + 0.5 * last_row['TIME INTERVAL'][0] return int(math.ceil(end_time)) ################################################################################ def check_table_version_exists(aips_data, table, version): """Check whether a table version exists, as 'AIPS CL', 'AIPS BP', etc. if version == 0, then check for ANY of that table, """ # If version is zero, then accept any of that table if(version == 0): for t in aips_data.tables: if(t[1] == table): break else: return 0 return 1 # Otherwise, check for the version number too for t in aips_data.tables: if((t[0] == version) and (t[1] == table)): break else: return 0 return 1 ################################################################################ def verify_table_version_free(aips_data, table, version, overwrite): """check that a table version is vacant, and possibly delete if ok""" # check for an existing table version if(check_table_version_exists(aips_data, table, version)): # Are we allowed to delete it? if(overwrite): aips_data.zap_table(table, version) else: raise RuntimeError, "Table %s %d exists in dataset '%s', but cannot delete"%(table, version, aips_data.__str__()) ################################################################################ def clear_table_version_maybe(aips_data, table, version, overwrite): """if we are overwriting and the table exists, zap it""" if(overwrite): # check for an existing table version if(check_table_version_exists(aips_data, table, version)): aips_data.zap_table(table, version) ################################################################################ # Get the Year, month, day for the start of observations def get_observation_year_month_day(aips_data): date_string = aips_data.header.date_obs date_list = date_string.split('-') year = int(date_list[0]) month = int(date_list[1]) day = int(date_list[2]) return (year, month, day) ################################################################################ # Get the reference day from the AN table def get_AN_reference_year_month_day(an_table): date_string = an_table.keywords['RDATE'] if(len(date_string) == 1): raise RuntimeError, "Error: The Obit string bug has finally been fixed" year = int(date_string[0:4]) month = int(date_string[4:6]) day = int(date_string[6:8]) return (year,month,day) ################################################################################ # Find the highest numbered antenna def find_antenna_highnumber(aips_data): """Find the highest numbered antenna in the dataset """ highnum = 0 for an in xrange(1,aips_data.table_highver('AIPS AN')+1): # check that this exists if not(check_table_version_exists(aips_data, 'AIPS AN', an)): continue # Ok, grab the AN table an_table = aips_data.table('AIPS AN', an) for row in an_table: if(highnum < row.nosta): highnum = row.nosta return highnum ################################################################################ # Find the highest numbered source def find_source_highnumber(aips_data): """Find the highest numbered source in the dataset """ highnum = 0 SU_table = aips_data.table('AIPS SU', 1) for row in SU_table: if(highnum < row.id__no): highnum = row.id__no return highnum ################################################################################ # Count visibilities def count_visibilities(aips_data_orig): """Count the number of visibilities per antenna in a dataset aips_data_orig I a non-Wizardry object to use """ aips_data = Wizardry.AIPSData.AIPSUVData(aips_data_orig.name, aips_data_orig.klass, aips_data_orig.disk, aips_data_orig.seq) # What is the highest valued antenna number? high_ant = find_antenna_highnumber(aips_data_orig) antennas = numarray.zeros((high_ant+1)) for row in aips_data: baseline = row.baseline a1 = baseline[0] a2 = baseline[1] if(a1 != a2): antennas[a1] = antennas[a1] +1 antennas[a2] = antennas[a2] +1 return antennas ################################################################################ # Get sources and positions def get_source_positions(aips_data): """Reads in source numbers, names, positions from SU table""" SU_table = aips_data.table('AIPS SU', 1) highnum = 0 for row in SU_table: if(highnum < row.id__no): highnum = row.id__no sources = (highnum+1)*[None] for row in SU_table: # verify that EPOCH is 2000.0 EPOCH = row.epoch if(EPOCH != 2000.0): raise RuntimeError, "Error: The AlbusIonosphere can only deal with J2000 coordinates right now." RA = row.raepo * M_DEG2RAD Dec = row.decepo * M_DEG2RAD name = row.source id = row.id__no sources[id] = [id,name,RA,Dec,EPOCH] # now return the table return sources ################################################################################ # set up the AlbusIonosphere for the AIPS reference information # and get antenna info def setup_ALbusIonosphere_for_ref_date(aips_data, AN_version): """set up the AlbusIonosphere for the AIPS reference information and get antenna info aips_data I the incoming UV dataset AN_version I the antenna table version to use. You probably want 1 0 should give the highest. OUTPUT antenna_info,subarray_info antenna_info O A list of antenna stuff, in mixed order, as [station_number,name,pos[x,y,z]] subarray_info O A list of keyword values for this subarray, as [year,month,day,hour, minute,second,POLARX,POLARY,UT1UTC] """ # Ok, grab the AN table an_table = aips_data.table('AIPS AN', AN_version) # Now, get the center of the array, as a numarray so we can do math # AIPS seems to be using a left-handed XYZ coordinate system, # SOMETIMES, for certain antennas. ARRRGHHH!!! array_name = an_table.keywords['ARRNAM'] if(len(array_name) == 1): raise RuntimeError, "Error: the keyword string bug is finally fixed" array_center = numarray.array([an_table.keywords['ARRAYX'], an_table.keywords['ARRAYY'], an_table.keywords['ARRAYZ']]) MAGIC_SIZE = 1000. MAGIC_MULTIPLIER = 1.0 if( ( (abs(array_center[0]) < MAGIC_SIZE) and (abs(array_center[1]) < MAGIC_SIZE) and (abs(array_center[2]) < MAGIC_SIZE) ) and not( (array_name == 'ATCA ') or (array_name == 'ATLBA ') ) ): # using left-handed system MAGIC_MULTIPLIER = -1.0 array_center[1] = array_center[1] * MAGIC_MULTIPLIER # AIPS dataset reference date (year_ref,month_ref,day_ref) = get_observation_year_month_day(aips_data) # The AN reference date (year_an, month_an, day_an) = get_AN_reference_year_month_day(an_table) if((year_ref,month_ref,day_ref) != (year_an, month_an, day_an)): #raise RuntimeError, "The reference dates don't match!" warnings.warn("The reference dates don't match!") # Now, get the time offset from the AIPS time to UTC UTC_offset = -an_table.keywords['DATUTC'] # Now, set the AlbusIonosphere reference time retval = AlbusIonosphere.set_reference_time_AIPS(year_ref, month_ref, day_ref, 0, 0, UTC_offset, an_table.keywords['POLARX'], an_table.keywords['POLARY'], an_table.keywords['UT1UTC']) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.set_reference_time_AIPS gave %d"%retval # Make a list of these parameters to send back to the caller subarray_info = [year_ref, month_ref, day_ref, 0, 0, UTC_offset, an_table.keywords['POLARX'], an_table.keywords['POLARY'], an_table.keywords['UT1UTC'] ] # Now make a list of the antenna station numbers and positions. antenna_info = [] for row in an_table: station = row.nosta pos = numarray.array(row.stabxyz) pos[1] = pos[1] * MAGIC_MULTIPLIER pos = pos + array_center name = row.anname if(row.mntsta == 2): warnings.warn("Orbiting Station '%s'found --- setting ionosphere to 0"%name) pos = pos * 0.0 antenna_info.append([station,name,pos]) return antenna_info,subarray_info ################################################################################ # make a proper array of station positions def make_proper_antenna_postion_array(antenna_info): """convert a Python antenna_info list to a numarray of positions""" #First, what is the maximum station number? max_num = 0 for ant in antenna_info: if(ant[0] > max_num): max_num = ant[0] # Now allocate an array antenna_array = numarray.zeros((max_num+1,3),type='Float64') # Now fill things in for ant in antenna_info: antenna_array[ant[0],:] = ant[2] # that's it return antenna_array ################################################################################ # Setup AlbusIonosphere for source def setup_ALbusIonosphere_for_source(aips_data, processing_option, source): """Setup AlbusIonosphere for source aips_data I the main uv dataset processing_option I the text string of the processing option source I Source information [ID, NAME, RA, Dec, EPOCH] """ # First, tell Albus stuff where the source is retval = AlbusIonosphere.set_source_position(source[2], source[3]) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.set_source_position gave %d"%retval # do anything else that needs to be done here, perhaps # based on the processing option. For now, this is all return ################################################################################ # Setup AlbusIonosphere for source def setup_ALbusIonosphere_for_station(aips_data, processing_option, station): """Setup AlbusIonosphere for source aips_data I the main uv dataset processing_option I the text string of the processing option station I Station information [ID, NAME, [x,y,z]] The return value of this function depends on whether the station is a valid Earth-based station. 1 means Earth-based, 0 means space-based """ # First, check for a strange antenna location pos = station[2] if(numarray.innerproduct(pos,pos) == 0.0): # strange station return 0 # Next, tell Albus stuff where the station is retval = AlbusIonosphere.set_station_position(pos[0],pos[1],pos[2]) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.set_station_position gave %d"%retval # do anything else that needs to be done here, perhaps # based on the processing option. For now, this is all return 1 ################################################################################ # Add blank entry to the AI table def add_blank_entry_to_AI(AI_table, antenna_id, source_id, subarray, time, time_width): """Add a blank entry to the AI table, to let a spacecraft observation through AI_table I The AI table to append to antenna_id I station number source_id I which source subarray I which subarray time I time in seconds since reference date time_width I full width in seconds for which solution to be blanked """ new_row = Wizardry.AIPSData.AIPSTableRow(AI_table) new_row.time = time new_row.time_interval = time_width new_row.source_id = source_id new_row.antenna_no = antenna_id new_row.subarray = subarray new_row.az = 0.0 new_row.el = 0.0 new_row.ion_data = [0.0,0.0,1.0] AI_table.append(new_row) ################################################################################ # process a scan of data for AlbusIonosphere information def process_AlbusIonosphere_for_scan(AI_table, antenna_id, source_id, subarray, start_time, end_time): """process a scan of data for AlbusIonosphere information This function assumes that everything except for a specific scan has been set up in the AlbusIonosphere information area. This function will run through all of the times the AlbusIonosphere makes data, and appends the information tot he AI_table AI_table I The AI table to append to antenna_id I station number source_id I which source subarray I which subarray start_time I start of the scan in seconds since the reference day end_time I end of the scan in seconds since the reference day """ # Tell AlbusIonosphere about the times retval = AlbusIonosphere.set_scan_times(start_time, end_time) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.set_scan_times gave %d"%retval # Get a blank row row = Wizardry.AIPSData.AIPSTableRow(AI_table) # Now, how many times do we process Num_Data = AlbusIonosphere.get_Num_Ionospheric_Predictions() for i in xrange(Num_Data): # Calculate the Ionosphere retval,pred_time, time_width, El, Az, STEC, SRM, VTEC_factor = \ AlbusIonosphere.get_ionospheric_prediction(i) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.get_ionospheric_prediction gave %d"%retval row.time = pred_time * DAYS_PER_SECOND row.time_interval = time_width * DAYS_PER_SECOND row.source_id = source_id row.antenna_no = antenna_id row.subarray = subarray row.az = Az * M_RAD2DEG row.el = El * M_RAD2DEG row.ion_data = [STEC, SRM, VTEC_factor] #print row AI_table.append(row) ################################################################################ # Add information from a subarray to the AI table def add_subarray_to_AI(aips_data, processing_option, AI_table, NX_table, AN_version, sources): """Add information from a subarray to the AI table aips_data I the main uv dataset processing_option I the text string of the processing option AI_table B The (opened) AI table to add information to NX_table I The (opened) NX table to get scan info from AN_version I The subarray version number to process sources I List of sources to get Ionosphere data for """ # First, get the antenna information for subarray AN_version antenna_info,subarray_info = \ setup_ALbusIonosphere_for_ref_date(aips_data, AN_version) # Loop over the antennas #print antenna_info for antenna in antenna_info: station_id = antenna[0] print "Working on Antenna %2d ('%s')"%(station_id,antenna[1]) # Cehck if we need to set up any special antenna-based information # in the AlbusIonosphere area station_code = setup_ALbusIonosphere_for_station(aips_data, processing_option, antenna) # Now, work our way through the scans for scan in NX_table: # If this is not our subarray, then skip if(scan.subarray != AN_version): continue source_id = scan.source_id # Get the start and stop times mid_time = scan.time * SECONDS_PER_DAY Delta_time = scan.time_interval * SECONDS_PER_DAY start_time = mid_time - 0.5 * Delta_time end_time = mid_time + 0.5 * Delta_time # If we are working on a normal station if(station_code): # Now do each source. If source_id is out of range, this # will throw an exception, which is what I want. source = sources[source_id] if(source != None): print " Working on Source %2d ('%s')"%(source_id,source[1]) # Check if we need to set up something special in ALbus setup_ALbusIonosphere_for_source(aips_data, processing_option, source) process_AlbusIonosphere_for_scan(AI_table, station_id, source_id, AN_version, start_time, end_time) else: # No data for this source raise RuntimeError, "Error: Cannot find source %d"%source_id else: add_blank_entry_to_AI(AI_table, station_id, source_id, AN_version,start_time-1.0,Delta_time+2.0) add_blank_entry_to_AI(AI_table, station_id, source_id, AN_version,end_time+1.0,Delta_time+2.0) # end of station_code check # for scan over NX table # for antenna over antenna info return ################################################################################ # Read in the data from an AI table def read_in_AI(aips_data, AI_version): """Read in the entire AI table to memory aips_data I an AIPS uv dataset (standard, not Wizardry) AI_version I version number of the AI table to create This function will return AI_header_data,AI_table_data where AI_header_data O the important parts of the table information NO_ANT, NO_SRC, NO_TERM, IONOTYPE, DATA_SRC AI_table_data O A list of the table rows """ # Ok, open the table AI_table = aips_data.table('AIPS AI', AI_version) num_rows = len(AI_table) AI_dataset = Wizardry.AIPSData.AIPSUVData(aips_data.name, aips_data.klass, aips_data.disk, aips_data.seq) AI_table = AI_dataset.table('AIPS AI', AI_version) AI_header_data = {} AI_header_list = ['NO_ANT', 'NO_SRC', 'NO_TERM', 'IONOTYPE', 'DATA_SRC'] for item in AI_header_list: AI_header_data[item] = AI_table.keywords[item] # Now read in all of the rows. Predeclare some numarray arrays to store # things in #AI_table_data = [] AI_table_data = {} time = numarray.zeros((num_rows), type='Float64') time_interval = numarray.zeros((num_rows), type='Float64') source_id = numarray.zeros((num_rows), type='Int32') antenna_no = numarray.zeros((num_rows), type='Int32') subarray = numarray.zeros((num_rows), type='Int32') az = numarray.zeros((num_rows), type='Float64') el = numarray.zeros((num_rows), type='Float64') ion_data = numarray.zeros((num_rows,AI_header_data['NO_TERM']), type='Float64') i=0 for row in AI_table: time[i] = row.time time_interval[i] = row.time_interval source_id[i] = row.source_id antenna_no[i] = row.antenna_no subarray[i] = row.subarray az[i] = row.az el[i] = row.el ion_data[i,:] = row.ion_data i = i+1 AI_table_data['time'] = time AI_table_data['time_interval'] = time_interval AI_table_data['source_id'] = source_id AI_table_data['antenna_no'] = antenna_no AI_table_data['subarray'] = subarray AI_table_data['az'] = az AI_table_data['el'] = el AI_table_data['ion_data'] = ion_data return AI_header_data,AI_table_data ################################################################################ # Find the upper and lower indices for interpolating in the AI table def find_AI_interpolation_high_low(AI_table, subarray_req, antenna_req, source_req, time_req, low, high): """Find the upper and lower indices for interpolating in the AI table AI_table I The dict of columns for the AI table. Note that this is NOT the AIPS table, but the dict of all columns. The entries should be sorted by time, source, antenna, so that all entries for a single source are together in increasing time order for a single telescope/subarray. See read_in_AI for the description of the table subrarray_req I The subarray value antenna_req I The antenna source_req I The source ID time _req I The time slot low I The starting guess for the low index 0 to initialize high I The starting guess for the high index Note that if low or high are not simple integers, they might get modified. """ assert(low >= 0) time = AI_table['time'] subarray = AI_table['subarray'] antenna = AI_table['antenna_no'] source = AI_table['source_id'] size = source.size() if(low == 0): while(low < size): if( (source[low] != source_req) or (antenna[low] != antenna_req) or (subarray[low] != subarray_req) ): low = low + 1 else: break # This is the first time for this source. Check that it is valid if((time[low] > time_req) or (low >= size)): # Not valid low = -1 high = -1 return low,high if(time[low] > time_req): # How did this slot have a time too great? Try again low = 0 return find_AI_interpolation_high_low(AI_table, subrarray_req, antenna_req, source_req, time_req, low, high) # Now, working from the last valid low, find others. low_last = low low = low + 1 while(low < size): if( (source[low] != source_req) or (antenna[low] != antenna_req) or (subarray[low] != subarray_req) ): # Not a valid antenna/source, so continue pass else: # This is a valid antenna/soure. Is thie time ok? if(time[low] <= time_req): # This is still a valid option low_last = low low = low + 1 else: # the time is too large, so stop break low = low + 1 low = low_last # Now do the high part if(high <= low): high = low+1 while(high < size): if( (time[high] <= time_req) or (source[high] != source_req) or (antenna[high] != antenna_req) or (subarray[high] != subarray_req) ): high = high + 1 else: break if(high >= size): # If we have ended here, then there is no data low = -1 high = -1 # low and high found return low,high ################################################################################ # Find the upper and lower indices for interpolating in the AI table def find_AI_nearest_time_point(AI_table, subarray_req, antenna_req, source_req, time_req, last_index): """Find the upper and lower indices for interpolating in the AI table AI_table I The dict of columns for the AI table. Note that this is NOT the AIPS table, but the dict of all columns. The entries should be sorted by time, source, antenna, so that all entries for a single source are together in increasing time order for a single telescope/subarray. See read_in_AI for the description of the table subrarray_req I The subarray value antenna_req I The antenna source_req I The source ID time _req I The time slot last_index I The last index as to where the nearest time point was found. Note that if this is not a simple integer, it might get modified. """ assert(last_index >= 0) time = AI_table['time'] subarray = AI_table['subarray'] antenna = AI_table['antenna_no'] source = AI_table['source_id'] size = source.size() # make sure that our index points to a valid antenna/source index = last_index while(index < size): if( (source[index] != source_req) or (antenna[index] != antenna_req) or (subarray[index] != subarray_req) ): index = index + 1 else: break if(index >= size): # if we have made it all the way through, then something is strange, and # start over from 0 if(last_index != 0): index=0 return find_AI_nearest_time_point(AI_table, subarray_req, antenna_req, source_req, time_req, index) else: # there were no valid points! return -1 # Ok, this is a valid point. Get the time difference Delta_t_min = abs(time_req - time[index]) # Check for going down in index last_index = index index = index-1 while(index >= 0): if( (source[index] == source_req) and (antenna[index] == antenna_req) and (subarray[index] == subarray_req) ): Delta_t = abs(time_req - time[index]) if(Delta_t > Delta_t_min): break else: Delta_t_min = Delta_t last_index = index index = index-1 # Check for increasing index index = last_index+1 while(index < size): if( (source[index] == source_req) and (antenna[index] == antenna_req) and (subarray[index] == subarray_req) ): Delta_t = abs(time_req - time[index]) if(Delta_t > Delta_t_min): break else: Delta_t_min = Delta_t last_index = index index = index+1 # Ok, this is it return last_index ################################################################################ # linearly interpolate ion data in STEC space def linear_interpolate_STEC(ion_data, time_array, index1, index2, calc_time): """linearly interpolate ion data in STEC space ion_data I a numarray[*,*] which holds the ionospheric data as STEC = ion_data[row,0] SRM = ion_data[row,1] VTEC_factor = ion_data[row,2] time_array I a numarray[*] which holds the ionosphere times as time_array[row] index1 I The lower time index index2 I The upper time index calc_time I The desired time OUTPUT STEC,SRM STEC O The slant total electron content, in electrons m^{-2} SRM 0 The slant rotation measure, in electrons m^{-2} T """ t1 = time_array[index1] t2 = time_array[index2] if(t1 == t2): # Hey, the times are the same! if(t1 == calc_time): # But they are both the needed time, so don't worry return ion_data[index1,0],ion_data[index1,1] else: raise RuntimeError, "Error: Both times are the same!" # Ok, calculate the slope factor factor = (calc_time - t1) / (t2-t1) STEC = ion_data[index1,0] + (ion_data[index2,0] - ion_data[index1,0]) *factor SRM = ion_data[index1,1] + (ion_data[index2,1] - ion_data[index1,1]) *factor #print STEC, SRM,t1,t2,calc_time #print ion_data[index1,0], ion_data[index2,0] return STEC,SRM ################################################################################ # linearly interpolate ion data in VTEC space def linear_interpolate_VTEC(ion_data, time_array, index1, index2, calc_time): """linearly interpolate ion data in STEC space ion_data I a numarray[*,*] which holds the ionospheric data as STEC = ion_data[row,0] SRM = ion_data[row,1] VTEC_factor = ion_data[row,2] time_array I a numarray[*] which holds the ionosphere times as time_array[row] index1 I The lower time index index2 I The upper time index calc_time I The desired time OUTPUT STEC,SRM STEC O The slant total electron content, in electrons m^{-2} SRM 0 The slant rotation measure, in electrons m^{-2} T """ t1 = time_array[index1] t2 = time_array[index2] if(t1 == t2): # Hey, the times are the same! if(t1 == calc_time): # But they are both the needed time, so don't worry return ion_data[index1,0],ion_data[index1,1] else: raise RuntimeError, "Error: Both times are the same!" # Ok, calculate the slope factor factor = (calc_time - t1) / (t2-t1) VTEC = ion_data[index1,0]*ion_data[index1,2] \ + (ion_data[index2,0]*ion_data[index2,2] - ion_data[index1,0]*ion_data[index1,2]) *factor VRM = ion_data[index1,1]*ion_data[index1,2] \ + (ion_data[index2,1]*ion_data[index2,2] - ion_data[index1,1]*ion_data[index1,2]) *factor V = ion_data[index1,2] + (ion_data[index2,2] - ion_data[index1,2]) *factor STEC = VTEC / V SRM = VRM / V #print V, factor, VTEC, VRM, STEC, SRM,t1,t2,calc_time #print ion_data[index1,0], ion_data[index2,0], ion_data[index1,2],ion_data[index2,2] return STEC,SRM ################################################################################ # linearly interpolate ion data in VTEC space using the actual source VTEC factor def linear_interpolate_VTEC_source(ion_data, time_array, index1, index2, calc_time, subarray_info_list, subarray, antenna_info_array, antenna, source_info_array, source): """linearly interpolate ion data in STEC space using the actual source VTEC factor ion_data I a numarray[*,*] which holds the ionospheric data as STEC = ion_data[row,0] SRM = ion_data[row,1] VTEC_factor = ion_data[row,2] time_array I a numarray[*] which holds the ionosphere times as time_array[row] index1 I The lower time index index2 I The upper time index calc_time I The desired time subarray_info_list I LIst of information about subarray stuff subarray I The current subarray to use antenna_info_list I list of numarray array of antenna positions antenna_info_list[subarray][antenna][XYZ] antenna I The current antenna number source_info_list I list of source position lists source I The current source number OUTPUT STEC,SRM STEC O The slant total electron content, in electrons m^{-2} SRM 0 The slant rotation measure, in electrons m^{-2} T """ # How much of the AlbusIonosphere stuff do we need to update? if(subarray != linear_interpolate_VTEC_source.last_subarray): # This means that everything must be updated linear_interpolate_VTEC_source.last_subarray = subarray linear_interpolate_VTEC_source.last_antenna = -1 linear_interpolate_VTEC_source.last_source = -1 s_list = subarray_info_list[subarray] #print 'Setting up for subarray', subarray, s_list retval = AlbusIonosphere.set_reference_time_AIPS(s_list[0], s_list[1], s_list[2], s_list[3], s_list[4], s_list[5], s_list[6], s_list[7], s_list[8]) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.set_reference_time_AIPS gave %d"%retval if(antenna != linear_interpolate_VTEC_source.last_antenna): linear_interpolate_VTEC_source.last_antenna = antenna pos = antenna_info_array[subarray] #print 'Setting up for antenna', antenna, pos[antenna,0], pos[antenna,1], pos[antenna,2] retval = AlbusIonosphere.set_station_position(pos[antenna,0], pos[antenna,1], pos[antenna,2]) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.set_station_position gave %d"%retval if(source != linear_interpolate_VTEC_source.last_source): # Get the soruce info. If this is out of range, then # this raises an exception. If it is None, raise my own source_info = source_info_array[source] if(source_info == None): raise RuntimeError, "Error: no data for source %d"%source linear_interpolate_VTEC_source.last_source = source #print 'setting up for source', source, source_info retval = AlbusIonosphere.set_source_position(source_info[2], source_info[3]) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.set_source_position gave %d"%retval # Now start dealing with the interpolation t1 = time_array[index1] t2 = time_array[index2] if(t1 == t2): # Hey, the times are the same! if(t1 == calc_time): # But they are both the needed time, so don't worry return ion_data[index1,0],ion_data[index1,1] else: raise RuntimeError, "Error: Both times are the same!" # Get the proper VTEC_factor retval, El, Az, VTEC_factor = \ AlbusIonosphere.get_source_AzElVTEC(calc_time * SECONDS_PER_DAY, 300.0E3) if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere.get_source_AzElVTEC gave %d"%retval # Ok, calculate the slope factor factor = (calc_time - t1) / (t2-t1) VTEC = ion_data[index1,0]*ion_data[index1,2] \ + (ion_data[index2,0]*ion_data[index2,2] - ion_data[index1,0]*ion_data[index1,2]) *factor VRM = ion_data[index1,1]*ion_data[index1,2] \ + (ion_data[index2,1]*ion_data[index2,2] - ion_data[index1,1]*ion_data[index1,2]) *factor STEC = VTEC / VTEC_factor SRM = VRM / VTEC_factor #print VTEC_factor, factor, VTEC, VRM, STEC, SRM,t1,t2,calc_time #print ion_data[index1,0], ion_data[index2,0], ion_data[index1,2],ion_data[index2,2] return STEC,SRM linear_interpolate_VTEC_source.last_subarray = -1 linear_interpolate_VTEC_source.last_antenna = -1 linear_interpolate_VTEC_source.last_source = -1 ################################################################################ # Get some info lists def get_subarray_antenna_source_info(aips_data): """Get some info lists Note that this function call will leave the AlbusIonosphere stuff in a pretty-much indetreminate state for station/reference date information. aips_data I A normal AIPSUVDAta object OUTPUT subarray_info O List of information about subarray stuff, each element is [year,month,day,hour, minute,second,POLARX,POLARY,UT1UTC] antenna_info O list of numarray array of antenna positions antenna_info_list[subarray][antenna][XYZ] from center of Earth, in a right-handed coordinate system source_info O list of source position lists l = source_info[source_number] gives a list l[0] = source_id number l[1] = name as a string l[2] = RA (radians) l[3] = Dec (radians) l[4] = EPOCH (years) """ # Create the top-level subarray and antenna lists num_subarrays = aips_data.table_highver('AIPS AN') subarray_info = (num_subarrays+1)*[None] antenna_info = copy.copy(subarray_info) # Ok, for all subarrays, add stuff on for an in xrange(1,num_subarrays+1): # check that this exists if not(check_table_version_exists(aips_data, 'AIPS AN', an)): continue # Ok, fill in the list stuff ant,sub = setup_ALbusIonosphere_for_ref_date(aips_data, an) ant_numarray = make_proper_antenna_postion_array(ant) subarray_info[an] = sub antenna_info[an] = ant_numarray source_info = get_source_positions(aips_data) return subarray_info,antenna_info,source_info ################################################################################ # Find the upper and lower indices for interpolating in the AI table def apply_AI_table_to_CL_table(aips_data, AI_version, gainver, gainuse, interpolation_method, overwrite=0): """Find the upper and lower indices for interpolating in the AI table aips_data I an AIPS uv dataset (standard, not Wizardry) AI_version I The version of the AI table to apply gainver I The CL table to use for input. Must exist gainuse I The CL table version to write out. Note that if overwrite is true, then this one will be deleted. If overwrite is false, then this table must not exist yet. interpolation_method I The method to use for interpolating between AI_table points to get CL table data. Valid options are 'NearestNeighbor' -- nearest point in time for the same antenna/source 'LinearInterpolateSTEC' -- linearly interpolate points in STEC space 'LinearInterpolateVTEC' -- linearly interpolate points in VTEC space 'LinearInterpolateVTECSource' -- linearly interpolate points in VTEC space using the proper VTEC_factor of the source at the time overwrite I flag to indicate whether the new CL table may be overwritten. """ # Typically, if gainuse is 0 then get the *next* highest if(gainuse == 0): gianuse = aips_data.table_highver('AIPS CL') + 1 # Next, make sure the new CL table is free verify_table_version_free(aips_data, 'AIPS CL', gainuse, overwrite) # Now read in the AI table print 'Reading in AI table', AI_version print 'This may take a long time. Try reading some e-mail.' start_time = systime.clock() AI_header,AI_table = read_in_AI(aips_data, AI_version) print 'Ok, I have read the AI table in %.1f seconds. Continuing.'%(systime.clock()-start_time) # Get some extra information if we really need it if(interpolation_method == 'LinearInterpolateVTECSource'): subarray_info,antenna_info,source_info = \ get_subarray_antenna_source_info(aips_data) # Get a Wizardy object W_dataset = Wizardry.AIPSData.AIPSUVData(aips_data.name, aips_data.klass, aips_data.disk, aips_data.seq) # Get the old table and the new one CL_table_in = W_dataset.table('CL', gainver) CL_table_out = W_dataset.attach_table( 'CL', gainuse, no_term=CL_table_in.keywords['NO_TERM']) # How big is the CL table? Num_Total_CL_Rows = len(CL_table_in) #no_if=CL_table_in.keywords['NO_IF'] This is automatically copied #no_pol=CL_table_in.keywords['NO_POL'] This is automatically copied CL_table_out.keywords['REVISION'] = CL_table_in.keywords['REVISION'] CL_table_out.keywords['NO_ANT'] = CL_table_in.keywords['NO_ANT'] CL_table_out.keywords['MGMOD'] = CL_table_in.keywords['MGMOD'] # how many polarizations are there? If 2, then set a flag two_pol = 0 if(CL_table_in.keywords['NO_POL'] > 1): two_pol = 1 # Create a couple of placeholders to store where we are in the AI_table num_ant = find_antenna_highnumber(aips_data) low = numarray.zeros((aips_data.table_highver('AIPS AN')+1,num_ant+1), type='Int32') high = numarray.zeros((aips_data.table_highver('AIPS AN')+1,num_ant+1), type='Int32') # directly hold the ion data and time information ion_data = AI_table['ion_data'] time_array = AI_table['time'] # Now loop over all of the rows in the old table to update and # place into the new one. Keep track of how far along I am, so that # I can print out some progress messages start_time = systime.clock() count = -1 for row in CL_table_in: # For testing, only do first few count = count + 1 #if(count > 10): # warnings.warn("Ending CL table loop eary for debugging") # break if( (count&0x3FF) == 0 ): print "%5.1f%% completed in %10.1f s"%\ ((float(count)*100.0/Num_Total_CL_Rows), systime.clock()-start_time) #if(count > 10000): # warnings.warn("Ending CL table loop eary for debugging") # break # Need information about this row subarray = row.subarray antenna = row.antenna_no source = row.source_id time = row.time STEC = 0.0 SRM = 0.0 # Which interpolation method to use? if(interpolation_method == 'NearestNeighbor'): index = find_AI_nearest_time_point(AI_table, subarray, antenna, source, time, low[subarray,antenna]) if(index >= 0): STEC = ion_data[index,0] SRM = ion_data[index,1] low[subarray,antenna] = index else: # No data, so apply a zero correction, which is already set # above low[subarray,antenna] = 0 elif(interpolation_method == 'LinearInterpolateSTEC'): index1,index2 = find_AI_interpolation_high_low(AI_table, subarray, antenna, source, time, low[subarray,antenna], high[subarray,antenna]) if(index1 >= 0): STEC,SRM = linear_interpolate_STEC(ion_data, time_array, index1, index2, time) low[subarray,antenna] = index1 high[subarray,antenna] = index2 else: # No data, so reset the positions low[subarray,antenna] = 0 high[subarray,antenna] = 0 elif(interpolation_method == 'LinearInterpolateVTEC'): index1,index2 = find_AI_interpolation_high_low(AI_table, subarray, antenna, source, time, low[subarray,antenna], high[subarray,antenna]) if(index1 >= 0): STEC,SRM = linear_interpolate_VTEC(ion_data, time_array, index1, index2, time) low[subarray,antenna] = index1 high[subarray,antenna] = index2 else: # No data, so reset the positions low[subarray,antenna] = 0 high[subarray,antenna] = 0 elif(interpolation_method == 'LinearInterpolateVTECSource'): index1,index2 = find_AI_interpolation_high_low(AI_table, subarray, antenna, source, time, low[subarray,antenna], high[subarray,antenna]) if(index1 >= 0): STEC,SRM = linear_interpolate_VTEC_source(ion_data, time_array, index1, index2, time, subarray_info,subarray, antenna_info, antenna, source_info, source) low[subarray,antenna] = index1 high[subarray,antenna] = index2 else: # No data, so reset the positions low[subarray,antenna] = 0 high[subarray,antenna] = 0 else: raise RuntimeError, "Error: unknown interpolation method '%s'"%interpolation_method #update the row disp = AIPS_TEC_TO_DISP * STEC rm = AIPS_RM_TO_RM * SRM row.i_far_rot = row.i_far_rot + rm row.disp_1 = row.disp_1 + disp if(two_pol): row.disp_2 = row.disp_2 + disp #print 'Final Values', disp, rm, STEC, SRM #print row CL_table_out.append(row) # Close the file to make sure the rows get written. CL_table_out.close() return ################################################################################ # negate the dispersive delay value def negate_disp_in_CL_table(aips_data, gainver, gainuse, overwrite=0): """negate the dispersive delay value aips_data I an AIPS uv dataset (standard, not Wizardry) gainver I The CL table to use for input. Must exist gainuse I The CL table version to write out. Note that if overwrite is true, then this one will be deleted. If overwrite is false, then this table must not exist yet. overwrite I flag to indicate whether the new CL table may be overwritten. """ # Typically, if gainuse is 0 then get the *next* highest if(gainuse == 0): gianuse = aips_data.table_highver('AIPS CL') + 1 # Next, make sure the new CL table is free verify_table_version_free(aips_data, 'AIPS CL', gainuse, overwrite) # Get a Wizardy object W_dataset = Wizardry.AIPSData.AIPSUVData(aips_data.name, aips_data.klass, aips_data.disk, aips_data.seq) # Get the old table and the new one CL_table_in = W_dataset.table('CL', gainver) CL_table_out = W_dataset.attach_table( 'CL', gainuse, no_term=CL_table_in.keywords['NO_TERM']) #no_if=CL_table_in.keywords['NO_IF'] This is automatically copied #no_pol=CL_table_in.keywords['NO_POL'] This is automatically copied CL_table_out.keywords['REVISION'] = CL_table_in.keywords['REVISION'] CL_table_out.keywords['NO_ANT'] = CL_table_in.keywords['NO_ANT'] CL_table_out.keywords['MGMOD'] = CL_table_in.keywords['MGMOD'] # how many polarizations are there? If 2, then set a flag two_pol = 0 if(CL_table_in.keywords['NO_POL'] > 1): two_pol = 1 # Now loop over the rows in the table if(two_pol): for row in CL_table_in: #update the row row.disp_1 = -row.disp_1 row.disp_2 = -row.disp_2 CL_table_out.append(row) else: for row in CL_table_in: #update the row row.disp_1 = -row.disp_1 CL_table_out.append(row) # Close the file to make sure the rows get written. CL_table_out.close() return ################################################################################ # Run a simple TECOF def ionosphere_simple_iononex(aips_data, processing_option, aparm=[None,1,0]): group_name = '' if(processing_option == 'IO_COD'): group_name='cod' elif(processing_option == 'IO_COR'): group_name='cor' elif(processing_option == 'IO_ESA'): group_name='esa' elif(processing_option == 'IO_ESR'): group_name='esr' elif(processing_option == 'IO_IGR'): group_name='igr' elif(processing_option == 'IO_IGS'): group_name='igs' elif(processing_option == 'IO_JPL'): group_name='jpl' elif(processing_option == 'IO_JPR'): group_name='jpr' elif(processing_option == 'IO_UPC'): group_name='upc' elif(processing_option == 'IO_UPR'): group_name='upr' else: raise RuntimeError, "unknown ionosphere processing option " + processing_option._str__() year = aips_data.header.date_obs[2:4] doy = get_observation_day_of_year(aips_data) filename = dirname + '/ionex/' + group_name + 'g%3.3d0.%si'%(doy,year) # Now run TECOR tecor = AIPSTask('tecor', version = aips_version_std) tecor.indata = aips_data tecor.infile = filename tecor.nfiles = get_number_days_observations(aips_data) tecor.gainver = Gainuse_Global tecor.gainuse = Gainuse_Global + 1 if(aparm != None): if(aparm[0] == None): tecor.aparm = aparm else: tecor.aparm[1:] = aparm tecor() ################################################################################ # Run a model Ionosphere to make an AI table def ionosphere_process_model(aips_data, processing_option, tolerance, AI_version, overwrite=1): """Make an AI table from a model ionosphere, using the AlbusIonosphere stuff aips_data I an AIPS uv dataset (standard, not Wizardry) processing_option I a string defining which model to use tolerance I fractional uncertainty in total electron content to ask for. Note that this is the requested fractional accuracy, not the achieved. For L-band EVN observations, you should probably ask for no better than 1E-4, and probably 1E-5. For lower frequencies, you want even better precision, probably 1E-6. AI_version I version number of the AI table to create overwrite I flag as to whether the given table version may be overwritten """ # check for an existing table version verify_table_version_free(aips_data, 'AIPS AI', AI_version, overwrite) # Start up the Ionosphere model retval = 0 if(processing_option == 'MO_IRI'): retval = AlbusIonosphere.set_ionosphere_IRI(tolerance,tolerance*100.0) elif(processing_option == 'MO_PIM'): retval = AlbusIonosphere.set_ionosphere_PIM(tolerance,tolerance*100.0) else: raise RuntimeError, "unknown ionosphere processing option " + processing_option._str__() if(retval < 0): raise RuntimeError, "Error: AlbusIonosphere initilization for '%s' gave %d"%(processing_option,retval) # Now open the table for output AI_dataset = Wizardry.AIPSData.AIPSUVData(aips_data.name, aips_data.klass, aips_data.disk, aips_data.seq) CL_table = AI_dataset.table('CL', 0) AI_table = AI_dataset.attach_table('AI', AI_version, no_term=3) AI_table.keywords['NO_ANT'] = CL_table.keywords['NO_ANT'] AI_table.keywords['NO_SRC'] = find_antenna_highnumber(aips_data) AI_table.keywords['IONOTYPE'] = 'TEST1'.ljust(8)[0:8] AI_table.keywords['DATA_SRC'] = processing_option.ljust(8)[0:8] # I need the NX table NX_table = AI_dataset.table('NX', 0) # I need the sources sources = get_source_positions(aips_data) #print sources # Ok, for all subarrays, add stuff on for an in xrange(1,aips_data.table_highver('AIPS AN')+1): # check that this exists if not(check_table_version_exists(aips_data, 'AIPS AN', an)): continue # Ok, fill in the table add_subarray_to_AI(aips_data, processing_option, AI_table, NX_table, an, sources) AI_table.close() return ################################################################################ ################################################################################ ################################################################################ # uvdata = AIPSUVData('N05C2', 'UVDATA', 1, 1) # cltable = uvdata.table('CL', 0) # aitable = uvdata.attach_table('AI', 0, no_term=6) # aitable.keywords['NO_ANT'] = cltable.keywords['NO_ANT'] # row = AIPSTableRow(aitable) # for ...: # row.time = time # aitable.append(row)