parseltongue:apply_ai_table_to_cl_table
function: apply_AI_table_to_CL_table
################################################################################ # 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
parseltongue/apply_ai_table_to_cl_table.txt · Last modified: 2007/07/09 15:19 by 127.0.0.1