; MPROF (IDL) ; C. David Brown ; Washington University ; ; Revised for PDS archive ; November 1999 ; Judd D. Bowman ; Washington University ;; ;This procedure reads and displays single-orbit MOLA PEDR files. ;The user can zoom in on profiles in an interactive plot window ;and create text files of the data. Second-order shot coordinate ;correction and frame crossover residual offset are included. MOLA ;PEDR version 2.7 is required. This procedure is written to read ;single-orbit PEDR files only. Use the SPLITPEDR procedure to break ;apart multiple-orbit PEDR files before running this procedure. ; ;Inputs: ; ;ORBIT -- The orbit number of interest. The procedure will use ;this number to locate a PEDR file in the default directory or the ;directory specified by the DIRECTORY keyword. ; ;DIRECTORY -- (Optional) Set this keyword to the full path of the ;directory containing the single-orbit PEDR files. End the path with ;the appropriate character for your operating system. This keyword ;supersedes the default directory in the code. ;On a PC: 'C:\temp\pedrs\single\' ;On Unix: '/home/user/pedrs/single/' ; ;WIDTH -- (Optional) Set this keyword to the width of the plot window ;in pixels. ; ;HEIGHT -- (Optional) Set this keyword to the height of the plot window ;in pixels ; ;Example command line calls: ; ;IDL>mprof, 1578, directory='c:\temp\pedrs\single\', width=1000, height=500 ;IDL>mprof, 3018, directory='/home/user/pedrs/single/' ;IDL>mprof, 20 pro mprof, $ orbit, $ directory=dir, $ width=nx, $ height=ny ;Working directory if not keyword_set(dir) then begin ;********************************************************************* ; USER: You can modify the line below to define a default directory ; used for searching for pedr files. Otherwise, you can use the ; DIRECTORY keyword to specify a directory when calling the procedure ;********************************************************************* dir = '/path/input_directory/' endif ;------------------------------------------------------------------------- ; DECLARE VARIABLES - Define constants used in the procedure and establish ; some variables ;------------------------------------------------------------------------- ;Default plot window dimensions if not keyword_set(nx) then nx = 800 if not keyword_set(ny) then ny = 300 ;Mean planetary radius (km) RADIUS = double(3388) ;Number of kilometers per degree on Mars DTOK = !DPI * RADIUS / double(180) ;Filler values for bad shots NODATA1 = 1.0d17 NODATA2 = NODATA1 / 1.0d6 ;A nifty trick to find the byte order of the local machine ;0 = PC (Little endian) ;1 = UNIX (Big endian) ENDIAN = byte(256, 0) ;Length of the PEDR header len_header = long(7760) ;Length of one PEDR record len_record = long(776) ;PEDR record structure (SIS V.2.7 12/12/98) pedr_rec = { $ ; START BYTE frame_time_sec: 0L, $ ; 001 frame_time_microsec: 0L, $ ; 005 orbit_reference_num: 0L, $ ; 009 areo_latitude: 0L, $ ; 013 areo_longitude: 0L, $ ; 017 radial_distance: 0L, $ ; 021 frame_midpt_range: 0L, $ ; 025 shot_quality_flag: bytarr(4), $ ; 029 shot_quality_desc_flag: bytarr(16), $ ; 033 shot_planet_radius: lonarr(20), $ ; 049 midpt_planet_radius: 0L, $ ; 129 r_ascension: 0L, $ ; 133 declination: 0L, $ ; 137 twist: 0L, $ ; 141 rcvd_pulse_energy: intarr(20), $ ; 145 reflect_transmit: intarr(20), $ ; 185 first_channel_number: bytarr(20), $ ; 225 rtnd_pulse_width: intarr(20), $ ; 245 rcvd_pulse_width: intarr(20), $ ; 285 para_delta_latitude: 0L, $ ; 325 para_delta_longitude: 0L, $ ; 329 crossover_residual: 0L, $ ; 333 frame_mid_latlon: lonarr(2), $ ; 337 transmit_power: intarr(20), $ ; 345 shot_class_code: intarr(20), $ ; 385 backgrnd_noise_counts: lonarr(8), $ ; 425 range_delay: 0L, $ ; 457 range_width: 0L, $ ; 461 chan_thresh_setting: intarr(8), $ ; 465 rcvr_chan_mask: 0, $ ; 481 alg_word_min_hits: 0, $ ; 483 alg_word_hit_cnt: 0, $ ; 485 frame_counter: 0, $ ; 487 trigger_channel: 0, $ ; 489 frame_index: 0, $ ; 491 pkt_source_hdr: bytarr(8), $ ; 493 pkt_time_sec: 0L, $ ; 501 pkt_time_msec: 0, $ ; 505 packet_fine_time: 0, $ ; 507 sub_com_area: bytarr(28), $ ; 509 orbit_qual_flag: 0, $ ; 537 attitude_flag: 0, $ ; 539 frame_local_time: 0, $ ; 541 phase_angle: 0, $ ; 543 solar_angle: 0, $ ; 545 emission_angle: 0, $ ; 547 atm_opacity: 0L, $ ; 549 frame_midpt_time: 0D, $ ; 553 raw_rcvd_pulse_energy: bytarr(20), $ ; 561 raw_rcvd_pulse_width: bytarr(20), $ ; 581 delta_sc_lat: 0L, $ ; 601 delta_sc_lon: 0L, $ ; 605 delta_sc_rad: 0L, $ ; 609 areoid_radius: 0L, $ ; 613 off_nadir_angle: 0L, $ ; 617 bits_per_shot: bytarr(20), $ ; 621 delta_areoid: 0L, $ ; 641 mola_clock_rate: 0L, $ ; 645 mola_range: lonarr(20), $ ; 649 range_correction: intarr(20), $ ; 729 delta_latitude: 0L, $ ; 769 delta_longitude: 0L $ ; 773 } ;------------------------------------------------------------------------- ; PRELIMINARIES - Search the specified directory for the PEDR FILE based ; on the orbit number given by the user when calling the procedure ;------------------------------------------------------------------------- ;Return all files containing the specified orbit number files = findfile(dir + 'ap' + string(orbit, format='(i5.5)') $ + '*', count=count) ;At least one file was found if count ne 0 then begin ;Take the first file found infile = files[0] ;No files were found endif else begin ;Update the log window print, 'MPROF: No input file for orbit ' + strtrim(orbit, 2) + '.' ;Exit the procedure return endelse ;------------------------------------------------------------------------ ; READ PEDR FILE - Open the input file and import each record into ; latitude, longitude, and elevation (hit) arrays ;------------------------------------------------------------------------ ;Update the log window print, 'MPROF: Reading input file (' + infile + ').' ;For the PC if ENDIAN eq 0 then begin ;Open the input file openr, unit, infile, /get_lun, /binary, /noautomode ;For UNIX endif else begin ;Open the input file openr, unit, infile, /get_lun endelse ;Skip past the input file header to the data records point_lun, unit, len_header ;Calculate the number of records in the input file ;based on the size of the file file_stats = fstat(unit) n_records = (file_stats.size - len_header) / len_record n_shots = n_records * 20 ;Create lat, lon, and hit arrays based on number of ;records in the input file lon = dblarr(n_shots) lat = dblarr(n_shots) hit = dblarr(n_shots) ;An array to hold the hit values for all shots in a record subhit = dblarr(20) ;An array of doubles for processing lat and lon darray = (dindgen(20) - 9.5) / 20.0 ;Shot tracker start_shot = long(0) ;Loop through each record in the input file for loop = long(0), n_records - 1 do begin ;Read the record readu, unit, pedr_rec ;Swap byte order if running on a PC if ENDIAN eq 0 then begin pedr_rec.shot_planet_radius = swap_endian(pedr_rec.shot_planet_radius) pedr_rec.midpt_planet_radius = swap_endian(pedr_rec.midpt_planet_radius) pedr_rec.frame_mid_latlon = swap_endian(pedr_rec.frame_mid_latlon) pedr_rec.shot_class_code = swap_endian(pedr_rec.shot_class_code) pedr_rec.areoid_radius = swap_endian(pedr_rec.areoid_radius) pedr_rec.delta_areoid = swap_endian(pedr_rec.delta_areoid) pedr_rec.delta_latitude = swap_endian(pedr_rec.delta_latitude) pedr_rec.delta_longitude = swap_endian(pedr_rec.delta_longitude) pedr_rec.para_delta_latitude = swap_endian(pedr_rec.para_delta_latitude) pedr_rec.para_delta_longitude = swap_endian(pedr_rec.para_delta_longitude) pedr_rec.crossover_residual = swap_endian(pedr_rec.crossover_residual) pedr_rec.orbit_qual_flag = swap_endian(pedr_rec.orbit_qual_flag) endif ;Increment the shot counter start_shot = 20 * loop ;Subtract mid point radius from the shot radius drad = double(pedr_rec.shot_planet_radius $ - pedr_rec.midpt_planet_radius) / 1.0d5 ;Get the shots with class code equal to zero inodat = where(pedr_rec.shot_class_code eq 0, count) ;Set the drad value to zero for the shots with code zero if (count ne 0) then drad(inodat) = double(0) ;Calculate intermediate sublats for all shots sublat1 = darray * double(pedr_rec.delta_latitude) sublat2 = drad * double(pedr_rec.para_delta_latitude) ;Calculate intermediate sublons for all shots sublon1 = darray * double(pedr_rec.delta_longitude) sublon2 = drad * double(pedr_rec.para_delta_longitude) ;Fill latitude and longitude values for all shots lat(start_shot:start_shot + 19) = (double(pedr_rec.frame_mid_latlon(0)) $ + sublat1 + sublat2) / 1.0d6 lon(start_shot:start_shot + 19) = (double(pedr_rec.frame_mid_latlon(1)) $ + sublon1 + sublon2) / 1.0d6 ;Calculate the intermediate subhit values for all shots geo0 = double(pedr_rec.areoid_radius) subgeo = darray * double(pedr_rec.delta_areoid) xover = double(pedr_rec.crossover_residual) subhit = double(pedr_rec.shot_planet_radius) - (geo0 + subgeo) - xover ;Set the subhit value to NODATA1 for shots with code zero if (count ne 0) then subhit(inodat) = NODATA1 ;Add the subhit to the hit array hit(start_shot:start_shot + 19) = subhit / 1.0d5 ;Repeat for next record in the input file endfor ;Close the input file free_lun, unit ;Update the log window if necessary if (pedr_rec.orbit_qual_flag lt 2) then begin print, 'MPROF: Unsupported PEDR version! Proceeding anyway.' endif ;Update the log window print, 'MPROF: ' + strtrim(n_records, 2) + ' records in input file.' print, 'MPROF: ' + strtrim(n_shots, 2) + ' total shots imported.' ;Find the minimum latitude ;If the index_min does not equal 0 or n_shots - 1 ;then the south polar region has been crossed minlat = min(lat, index_min) ;Find maximum latitude ;If the index_max does not equal 0 or n_shots - 1 ;then the north polar region has been crossed maxlat = max(lat, index_max) ;The latitude values go in descending order if index_min gt index_max then begin ;Reverse the order for loop = long(0), (n_shots / 2) - 1 do begin ;Latitude temp = lat[loop] lat[loop] = lat[n_shots - 1 - loop] lat[n_shots - 1 - loop] = temp ;Longitude temp = lat[loop] lon[loop] = lon[n_shots - 1 - loop] lon[n_shots - 1 - loop] = temp ;Elevation (hit) temp = hit[loop] hit[loop] = hit[n_shots - 1 - loop] hit[n_shots - 1 - loop] = temp endfor ;Reverse the index pointers index_min = n_shots - 1 - index_min index_max = n_shots - 1 - index_max endif ;Adjust latitude values for south polar crossings for loop = long(0), index_min do begin lat[loop] = 2 * minlat - lat[loop] endfor ;Adjust latitude values for north polar crossings for loop = index_max, n_shots - 1 do begin lat[loop] = 2 * maxlat - lat[loop] endfor ;No south polar crossing if (index_min eq 0) or (index_min eq n_shots - 1) then minlat = -90.0 ;No north polar crossing if (index_max eq 0) or (index_max eq n_shots - 1) then maxlat = 90.0 ;Remember a representative spacing of the shots for showing ;data points in plots dlat = abs(lat[n_shots / 2] - lat[(n_shots / 2) - 1]) ;----------------------------------------------------------------------- ; CREATE WINDOW - Done importing the pedr file. Continue by setting up the ; plot window on the screen ;----------------------------------------------------------------------- ;Update the log window with instructions for using the plot print, '----------------------------------------------' print, 'MPROF: Instructions for using the plot window:' print, '----------------------------------------------' print, '* Click and drag with the LEFT mouse button at' print, 'any time to zoom in on a portion of the plot.' print, '' print, '* In the full profile plot, click the CENTER' print, 'mouse button to exit the procedure. Do not close' print, 'the window manually.' print, '' print, '* In a zoom view, click the CENTER mouse button' print, 'to create a postscript file of the plot and a' print, 'text file of the data displayed in the plot.' print, '' print, '* Use CONTROL + LEFT mouse button to simulate a' print, 'CENTER button if your mouse does not have one.' print, '' print, '* Click the RIGHT mouse button to return to the' print, 'plot of the full profile. Use the RIGHT mouse' print, 'button to restart the plot after resizing the' print, 'the plot window on the screen.' print, '----------------------------------------------' ;For the PC if ENDIAN eq 0 then begin ;Prepare the plot window set_plot, 'win' ;Set the desired font defont = 'HELVETICA*THIN*24' ;For UNIX endif else begin ;Prepare the plot window set_plot, 'x' ;Set the desired font defont = '-adobe-new century schoolbook' $ + '-bold-r-normal--18-180-75-75-p-113-iso8859-1' endelse ;Open the plot window on the screen window, /free, title='MPROF - MOLA Single Orbit Profile Plot', $ xsize=nx, ysize=ny, colors=256 win_unit = !d.window ;Set the device characteristics ;RETAIN=1: X-WINDOWS BACKING STORE (SUN) ;RETAIN=2: IDL BACKING STORE (SGI and PC) device, retain=2, decomposed=0, set_graphics = 3 ;Apply the desired font ;device, set_font=defont ;Setup the constant plot characteristics !order = 1 !x.thick = 1.0 !y.thick = 1.0 !p.thick = 1.0 !p.charthick = 2.0 ;Set the specific plot characteristics - ;These define the portion of the profile to be plotted index_start = long(0) index_stop = n_shots - 1 x_start = lat(index_start) x_stop = lat(index_stop) y_start = float(-10) y_stop = float(20) ;----------------------------------------------------------------------- ; USER EVENT LOOP. - Responding to user events ;----------------------------------------------------------------------- ;Setup variables for event processing ;The exit flag signals to stop processing events and ;exit the procedure. Acceptable values are 0 and 1. f_exit = 0 ;The draw flag is used to redraw the plot. Acceptable ;values are 0 and 1. f_draw = 1 ;The level flag keeps track of the number of zooms ;performed by the user. Acceptable values are ;positive integers. f_level = 0 ;Start the event loop while (f_exit eq 0) do begin ;--------------------------------------------- ; DRAW PLOT - Redraw the entire plot if the ; draw flag is set to 1 ;--------------------------------------------- if (f_draw eq 1) then begin ;Reset the draw flag f_draw = 0 ;Draw the plot without showing the data - ;This draws the title and Y-axis and establishes ;coordinates for the window in the data frame plot, [0], [0], $ /nodata, $ xtitle='Latitude (degrees)', $ xrange=[x_start, x_stop], $ xstyle=5, $ ytitle='Elevation (km)', $ yrange=[y_start, y_stop], $ ystyle=1, $ title=infile, $ max_value=NODATA2 ;Get the upper left corner of the plot region in device coordinates c_left = convert_coord(x_start, y_start, /data, /to_device) ;Get the minlat region in device coordinates c_minlat = convert_coord(minlat, y_start, /data, /to_device) ;Get the maxlat region in device coordinates c_maxlat = convert_coord(maxlat, y_start, /data, /to_device) ;Get the lower right corner of the plot region in device coordinates c_right = convert_coord(x_stop, y_stop, /data, /to_device) ;Set the default center plot boundaries to cover the entire plot region c_center_left = c_left[0] c_center_right = c_right[0] ;Set the default center data boundaries to cover the entire plot region d_center_left = x_start d_center_right = x_stop ;The south polar region is in the plot if (c_minlat[0] gt c_left[0]) and $ (c_minlat[0] lt c_right[0]) then begin ;Draw the left portion of the X-axis plot, [0], [0], $ /nodata, $ /noerase, $ xrange=[2 * minlat - x_start, minlat], $ xstyle=1, $ yrange=[y_start, y_stop], $ ystyle=5, $ position=[c_left[0], c_left[1], c_minlat[0], c_right[1]], $ /device ;Reset the data coordinates - does not show in plot plot, [0], [0], $ /nodata, $ /noerase, $ xrange=[x_start, x_stop], $ xstyle=5, $ yrange=[y_start, y_stop], $ ystyle=5, $ position=[c_left[0], c_left[1], c_right[0], c_right[1]], $ /device ;Draw the dashed line oplot, [minlat, minlat], !y.crange, $ linestyle=2 ;Move the center X-axis left boundary c_center_left = c_minlat[0] ;Move the center plot left boundary d_center_left = minlat endif ;The north polar region is in the plot if (c_maxlat[0] gt c_left[0]) and $ (c_maxlat[0] lt c_right[0]) then begin ;Draw the right portion of the X-axis plot, [0], [0], $ /nodata, $ /noerase, $ xrange=[maxlat, 2 * maxlat - x_stop], $ xstyle=1, $ yrange=[y_start, y_stop], $ ystyle=5, $ position=[c_maxlat[0], c_maxlat[1], c_right[0], c_right[1]], $ /device ;Reset the data coordinates - does not show in plot plot, [0], [0], $ /nodata, $ /noerase, $ xrange=[x_start, x_stop], $ xstyle=5, $ yrange=[y_start, y_stop], $ ystyle=5, $ position=[c_left[0], c_left[1], c_right[0], c_right[1]], $ /device ;Draw the dashed line oplot, [maxlat, maxlat], !y.crange, $ linestyle=2 ;Move the center X-axis right boundary c_center_right = c_maxlat[0] ;Move the center plot right boundary d_center_right = maxlat endif ;Draw the center portion of the X-axis plot, [0], [0], $ /nodata, $ /noerase, $ xrange=[d_center_left, d_center_right], $ xstyle=1, $ yrange=[y_start, y_stop], $ ystyle=5, $ position=[c_center_left, c_left[1], c_center_right, c_right[1]], $ /device ;Reset the data coordinates - does not show in plot plot, [0], [0], $ /nodata, $ /noerase, $ xrange=[x_start, x_stop], $ xstyle=5, $ yrange=[y_start, y_stop], $ ystyle=5, $ position=[c_left[0], c_left[1], c_right[0], c_right[1]], $ /device ;Draw the X-axis label xyouts, 0.50, 0.03, 'Latitude (degrees)', $ /norm, $ align=0.5 ;Determine if the data points should be marked or if ;only a line should be drawn - based on amount of zoom if (x_stop - x_start) gt (100 * dlat) then begin ;Only draw the line oplot, lat[index_start:index_stop], hit[index_start:index_stop], $ max_value=y_stop endif else begin ;Draw the line and mark the data points oplot, lat[index_start:index_stop], hit[index_start:index_stop], $ max_value=y_stop, $ psym=-6, $ symsize=0.4 endelse endif ;--------------------------------------------- ; WAIT FOR THE USER TO DO SOMETHING ;--------------------------------------------- ;Reset the system mouse-state variable !mouse.button = 0 ;Wait for a mouse button to be pressed cursor, x_o_stop, y_o_stop, /data ;The user clicked the LEFT mouse button ;Hande a zoom event if (!mouse.button eq 1) then begin ;Set the window to XOR draw mode device, set_graphics = 6 ;Set the start coordinates to the location of ;the mouse click x_start = x_o_stop y_start = y_o_stop ;Draw a box as the user moves the mouse. Repeat until ;the user releases the mouse button while (!mouse.button ne 0) do begin ;Wait for the user to move the mouse cursor, x_stop, y_stop, /data, /change ;Draw the new box plots, [x_start, x_start, x_stop, x_stop, x_start], $ [y_start, y_stop, y_stop, y_start, y_start], $ /data ;Erase the old box plots, [x_start, x_start, x_o_stop, x_o_stop, x_start], $ [y_start, y_o_stop, y_o_stop, y_start, y_start], $ /data ;Remember the location of the current box x_o_stop = x_stop y_o_stop = y_stop ;Update the box while the user drags the mouse endwhile ;Set the window back to regular draw mode device, set_graphics = 3 ;The box has valid dimensions if ((x_start ne x_stop) and (y_start ne y_stop)) then begin ;Set flag for updating plot f_draw = 1 ;Increase the level flag f_level = f_level + 1 ;Put the horizontal box coordinates in order if (x_start gt x_stop) then begin ;Swap the start and stop coordinates so ;the start coordinate is less than the stop ;coordinate x_o_stop = x_stop x_stop = x_start x_start = x_o_stop endif ;Put the vertical box coordinates in order if (y_start gt y_stop) then begin ;Swap the start and stop coordinates so ;the start coordinate is less than the stop ;coordinate y_o_stop = y_stop y_stop = y_start y_start = y_o_stop endif ;Find the chunk of the arrays to plot based on the new ;dimensions for the zoom ;Pointers into the arrays index_start = long(0) index_stop = long(0) ;Find the start of the plot while (lat[index_start] lt x_start) and $ (index_start lt (n_shots - 1)) do begin ;increment the index value index_start = index_start + 1 endwhile ;Start the search for the end of the plot index_stop = index_start while (lat[index_stop] lt x_stop) and $ (index_stop lt (n_shots - 1)) do begin ;increment the index value index_stop = index_stop + 1 endwhile ;The box has invalid dimensions endif else begin ;Set flag to skip updating plot f_draw = 0 endelse ;The user clicked the MIDDLE mouse button ;Handle an exit event or output event endif else if (!mouse.button eq 2) then begin ;Wait until the user releases the mouse button cursor, x_o_stop, y_o_stop, 4 ;The zoom level is 0 so exit the procedure if (f_level eq 0) then begin ;Set the exit flag f_exit = 1 ;The zoom level is greater than 0 so output a text file ;of the data in the plot currently on the screen endif else begin ;Ge the output file name from the user outfile = dialog_pickfile(path=dir, /write) ;Make sure a file was specified if outfile ne '' then begin ;Find the valid data points valid = where(hit[index_start:index_stop] lt NODATA2, count) ;Only write the file if there are valid data points if (count gt 0) then begin ;Open the output file openw, out_unit, outfile, /get_lun ;Write the output header printf, out_unit, 'Created by MPROF on ' + systime(0) printf, out_unit, 'MOLA data from orbit: ' + strtrim(orbit, 2) printf, out_unit, 'Number of shots in file: ' + strtrim(count, 2) printf, out_unit, '' printf, out_unit, 'Longitude.E', 'Latitude.N', 'Elevation (m)', $ format = '(3A20)' ;Write the output data to the file for i = long(0), count - 1 do begin ;Convert the latitudes back from the ajustments made ;for plotting around polar regions ;Convert for wrap around south polar region if (lat[index_start + valid[i]] lt minlat) then begin latout = 2 * minlat - lat[index_start + valid[i]] ;Convert for wrap around notrh polar regin endif else if (lat[index_start + valid[i]] gt maxlat) then begin latout = 2 * maxlat - lat[index_start + valid[i]] ;Make no change endif else begin latout = lat[index_start + valid[i]] endelse ;Write the data to the file printf, out_unit, $ lon[index_start + valid[i]], $ latout, $ 1000 * hit[index_start + valid[i]], $ format = '(3F20.5)' ;Repeat for each valid shot in plot endfor ;Close output file free_lun, out_unit ;Update the log window print, 'MPROF: Output file (' + outfile $ + ') written successfully.' endif else begin ;Update the log window print, 'MPROF: Output file NOT created' $ + ' - no valid data in plot.' endelse endif else begin ;Update the log window print, 'MPROF: Output file NOT created' $ + ' - no file specified.' endelse endelse ;The user clicked the RIGHT mouse button ;Handle a redraw full profile plot event endif else if (!mouse.button eq 4) then begin ;Wait until the user releases the mouse button cursor, x_o_start, y_o_start, 4 ;Reset the level flag f_level = 0 ;Set the draw flag f_draw = 1 ;Reset the plot characteristics index_start = long(0) index_stop = n_shots - 1 x_start = lat[index_start] x_stop = lat[index_stop] y_start = float(-10) y_stop = float(20) endif ;Repeat loop endwhile ;------------------------------------------------------------------------- ; EXIT - Delete the window and exit the procedure ;------------------------------------------------------------------------- ;Delete the plot window wdelete, win_unit ;Update the log window print, 'MPROF: Done.' ;Exit the procedure end