; CENERGY.PRO ; ;Written By Aaron Ridley now at University of Michigan ;Testing, support and scientific advise by Juan Rodriguez, Phillips Lab ; (at Ball Aerospace Corp., Boulder, CO as of 1998) ;(Note: In Oct 1997, the name of A.F. Phillips Lab., Geophysics Directorate ; was changed to Air Force Research Lab., Space Vehicles Directorate, ; Battlespace Environment Division or AFRL/VSB for short) ; ; Completed in August 1994 ; Modified on January 6, 1995 - Jean H. James ; Orginal version only worked on VMS. Changed it to work on unix. ; Further modified on June 8, 1995 - Jean H. James ; Added ssm processing on March 5, 1996 ; Modified January-July 1998 - F. J. Rich ; Replaced FORTRAN subroutines to read data with IDL procedures ; Bug repair: ; Sept 28, 1998 - Fixed plot alignment error due to data gaps. - fjr ; March 8, 1999 - Fixed error in routines which read dm and ssj4 ; data which caused some data to be ignored (especially late in ; when there are large data gaps. - fjr ; Modified November 2002 - F. J. Rich ; Re-sized and re-positioned items so that they all fit within ; the plot screen (interactive mode) or within the bounding box ; (postscript mode). Added option for encapsulated postscript. ; Bug repair/Modification: Dec 2002 ; Code for plotting lines and labels had become non-functional. ; Fixed code to make this option work again. Also added the ; beginning of code for plotting J4/J5,DM,Ni(SM) and SSM data ; all on the same plot. ; Modification: Jan 2003 ; Code for plotting J4/J5,DM,Ni(SM) and SSM data on same plot is ; finished (plot option 4). ; Modification: Dec 2003 - February 2004 ; Code added for plotting SSJ5 data and removing background due ; to MeV particles which pass through the s/c skin and instrument ; case and get counted. The SSJ5 is approximately 100 times more ; sensitive to they particle. ; Modification: August-September 2004 ; Added new geometric factors which correct errors in previous ; table for F11-F15 and includes an aging factor for each of the ; four detectors in SSJ4. No aging factor for F16 yet. Add ability ; to use SSIES3 data in EDR file format. Added code to handle time ; code in old SSJ4 files ; Modification: January 2004 ; Bug fixes related to plotting gaps in SSM data and printing ; zeros when no ephemeris values are available ; Modification: March 2007 ; Changed SSIES3 logic to accept new (as of Oct 06) name of EDR file. ; Minor bug fix for multi-plot option: 28 November 2007 ; ; Please report bugs to: ; Frederick Rich, AFRL/VSBXP, Hanscom AFB, MA 01731 ; e-mail: frederick.rich@hanscom.af.mil ; ;------------------------------- ; Procedures with this program: ; ;pro particle, j4file$, begtime, x_length, imulf, eflux, ecnt, iflux, $ ; icnt, satel, mlat, mlt, j4op, j4cl, $ ; minimum_cglat, nsq, totlen, plotparsw, removebackground, iyear ;PRO create_raw2count,raw2count ;PRO read_geomfac, isatnumber, e_energy, e_intenergy, e_intnumber, e_difnumber, $ ; e_geomf, i_energy, i_intenergy, i_intnumber, i_difnumber, $ ; i_geomf, iyear ;PRO GETJGF,IDSAT,NYR,GME,GMI,EE,EI,IFLG,ZF,CF ;PRO raw2ephem,bdata,jday,hour,minute,isecond,iyear,glat,glon,alt,gglat110,gglon110,$ ; cglat110,mlthr,mltmin,mltsec ;PRO convert2counts, bdata, raw2count, counts, isatnumber, mode_flag, modeb_cycle,$ ; modeb_zone ;PRO convert2flux, counts, num_min_saved_for_plot, $ ; e_fluxes, i_fluxes, e_geomf, i_geomf, isatnumber ;PRO jbackgroundremoval,jcntin,jcntout,isat,abkg,itype ;PRO open_output_device, outtype, psfile,datafile_or_keyboard_input,$ ; ps_filename_problem ;PRO get_file_names, partfile$, removebackground, driffile$, ssmffile$, smfile$,$ ; psfile, datafile_or_keyboard_input, output_type, dm_or_ssm_option, $ ; inputerror, ps_encapsulated_option ;PRO get_time_int, dayofmon, begtime, endtime, length, satel,cor, ycol, yfil, $ ; maxvel, max_nt, minlat, nsq, nodm, nossm, vopt,$ ; datafile_or_keyboard_input, portrait_or_land_option, diswitch, flipi,$ ; vcount, vlabc, vstr, vlabx, vlines, intsw, imul, output_type ;PRO tick_stuff, begtime, length, nticks, tickval, ticklab, bt, et, imul ;PRO set_color_palette, output_type, ycol ;PRO plot_electron, ele2d, bt, et, nticks, tickval, ticklab, yticklab, change,$ ; plotparsw, dm_or_ssm_option, ycol ;PRO plot_ion, ion2d, bt, et, nticks, tickval, ticklab, yticklab, $ ; change, flipi, plotparsw, dm_or_ssm_option,ycol ;PRO plot_dm, dmdata, bt, et, nticks, tickval, ticklab, $ ; yb, ye, ytick, ytickval, yticklab, mlt, $ ; mlat, begtime, npts, cor, imul, nodm, ycol, $ ; change, vopt ;PRO segmentplot,tarray,yarray,ye,npts,lnstyle,icolor,samples_per_second ;PRO plot_ssm,ssmdata, bt, et, nticks, tickval, $ ; yb, ye, ytick, ytickval, yticklab, $ ; begtime, nossm, ycol, $ ; change, output_type, dm_or_ssm_option ;PRO PLOT_SMNI,smnidata, bt, et, nticks, tickval, ticklab, $ ; begtime, nosmni, change, ycol ;PRO plot_ee, eeflx, eave, ieflx, iave, ticklab, intsw, $ ; nticks, tickval, bt, et, yfil, yticklab, change ;PRO plot_color_bar, cyb, cye, ctick, ctickval, cticke, cticki, change, $ ; diswitch, output_type, dm_or_ssm_option,ycol ;PRO plot_extra, begtime, monname, satel, change, psfile, dm_or_ssm_option, $ ; bt, et, nticks, tickval, ticklab, imul, mlt, mlat, $ ; yb, ye, intsw, output_type, removebackground ;PRO add_points, vlines, pts, bt, et, vcount, pcount,$ ; vstr, vlabx, vlabc, change, baset, datafile_or_keyboard_input, $ ; dm_or_ssm_option ;PRO rave, avearr, nadd, stdev ;PRO interpo, eflx, iflx, npts, ele2d, ion2d, diswitch, output_type, dm_or_ssm_option ;PRO calc_flx, ecnt, icnt, isat, eeflx, ieflx, eave, iave, eflx, iflx, $ ; oeeflx, oieflx, enflx, inflx, deeflx, deiflx, output_type, iyear, $ ; begtime ;PRO dm_ies3, edrfilenm, cor, begtime, inter, dmdata, isatin, $ ; mlto, mlato, glato, glono, mlatsso, alto, minlat, nsq,$ ; nodmdata, dm_or_sm_datafound ;PRO rdapgafile,filenm,error_flag ;PRO dm, filenm$,cor, begtime, inter, dmdata, isatin, $ ; mlto, mlato, glato, glono, mlatsso, alto, minlat, nsq,$ ; nodmdata, dm_or_sm_datafound ;PRO decode_1min_dm, bdata, vvert, vhorz, shkp2, svap,isat,year,jday,hour,minute,$ ; geolat, geolong, maglat, mlt, maglon, glatsol, glonsol,$ ; glat110, glon110, mlat110, mlon110, invlat,alt1,alt2,$ ; bx,by,bz,ex,ey,ez,ssenpot,svbias,svip,srepel,sifree ;PRO COROTFIX, FLOW, nsec_saved, NMINSV, gmt, gglat, alt, ecos ;PRO ssm,smfilen$,begtime,inter,imulf,ssmdata, $ ; glat0, glon0,mlatss0,alt0,minlat, $ ; nsq, nossmdata,dm_or_sm_datafound, arrcglalo ;pro CGLALO, a, OLDLAT, OLDLON, CGLAT, CGLON ;pro sm_ni,smfile,begtime,inter,smnidata ;PRO decode_1min_sm, bdata, nsm, smpower, smdensity, smrms,isat,$ ; year,jday,hour,minute, geolat, geolong, maglat,$ ; mlt, maglon, glatsol, glonsol, glat110, glon110, mlat110, $ ; mlon110, invlat,alt1,alt2, bx,by,bz,ex,ey,ez,ssenpot,$ ; svbias,svip,srepel,sifree ;pro sm_ni_ies3,edrfilenm,begtime,inter,smnidata ; ---- MAIN PROGRAM ----- ;--------------------------- ; ; This is an IDL Program which reads DMSP Ion Drift Meter (IDM) data from ; the SSIES instrument and precipitating electron and ion data from the ; SSJ4 or SSJ5 instrument, then plots them to a postscript or HP print file or ; interactively to a MS-Window, X-Windows. ; ; ; The data files are unformatted binary files, which IDL can read. For ; VMS, the record length must match the record length in the IDL 'open' ; statement. ; ; The code accepts inputs to run in the following ways: ; ; 1. Plot one set of data on one page. ; ; Each plot can have up to 30 minutes on it. If the ending time is ; 30 minutes or less past the beginning time and the maximum time ; span for each plot is set for that difference, and the minimum ; latitude is set for 0.0, the program will plot one data set on ; one page. ; ; 2. Plot a number of plots in a row (consecutive times). ; ; Set the beginning time for the first time, set the ending time ; for the end of the last plot, then decide the maximum amount ; of time for each plot. The program will read in one set of data, ; plot it out, then read in another set and plot, ect., until the ; last time is equal or greater then the ending time. The minimum ; latitude should be set to 0.0 also. ; ; 3. Plotting a number of plots above a certain latitude. ; ; Same as 2, but with a minimum latitude set. The program will ; ask whether which hemisphere to plot, or both. ; ; The multiple plots are kept track of in the following way: ; ; The user inputs BEGTIME, which contains the following information: ; BEGTIME(0) = year ; BEGTIME(1) = month ; BEGTIME(2) = day of month ; BEGTIME(3) = beginning hour ; BEGTIME(4) = beginning minute ; BEGTIME(5) = julian day ; ; ENDTIME is the time in seconds which the user has chosen to be the ; last time plotted. ; CURTIME is the beginning time of the current plot, in seconds. ; ; After the plot is made, CURTIME is incremented by the time length ; of the plot and compaired to ENDTIME. If CURTIME >= ENDTIME, the ; program stops. ; ; NEW FEATURES: ; ; -Each data gathering subroutine can be used to determine the minimum ; Latitude. ; ; -Also, the program is set up in a way in which if the driftmeter data ; is PRESENT, but does not fall within the invariant letitude selected, ; the program will skip the Particle gathering subroutine. If the ; driftmeter data is not present (as in missing orbit), the program ; will search through the particle data to try to find the "missing" ; data. ; ; -A seperate plot is made of the orbit in polar coordinate, so the path ; of the satellite can be traced. If driftmeter data is present, 1 second, ; or more (if the data is at time intervals greater than 1 second), ; averaged cross track velocities are plotted along the satellite's ; path. ; This subroutine uses invarient latitude and magnetic local time, using ; linear interpolation to determine the location of the satellite ; between ephemorous data points. There are two "tricky" points to this: ; namely over the equator and over the poles. In order for the plots to ; be as acurate as posible, two different routines are used for the ; linear interpolation. One method (the equator) simply linear interpolates ; the mlat and mlt. The other method (the pole) transforms the mlat and ; mlt to cartesian coordinates and interpolates. The reason for this ; is that mlt near the pole is not linear, but the cartesian coordinates ; are. Mlat over the equator, if transformed into cartesian coordinate, ; skips the near equator region, so using the actual values is the best ; method. ; ; -The option for an ascii output file was added. This file contains ; one second (or more than one second, as in 4 second, depending on the ; time interval between driftmeter data points. One second is the minimum) ; data. The file outputs time, glat, glon, alt, mlat, mlt, cross track ; velocity, vertical velocity, standard deviations for horizontal and ; vertical velocities, total energy for electron, average energy for ; electrons, total energy for ions, average energy for ions. The first ; grouping of particle data is for IONOSPHERIC INPUT (incedent), meaning ; that the energies less than .44 keV are ignored. The second grouping of ; particle data is the actual observations AT THE SATELLITE. The glat ; and glong are interpolated in the same way as the mlat and mlt. If the ; time in between driftmeter measurements is more than 3 second, the ; energies are average. If less than 3 seconds, the energy which ; corresponds to the time is output. ; ; -Users can enter lines and text on plots. The lines are from the bottom ; of the driftmeter data to the top of the energy flux plot, and are ; drawn at a specified time (hour, minute, second). The text is displayed ; below the driftmeter data and above the ephemeris data. It is centered ; on a specified time and is 70% the size of the normal text. ; ;----------------------------------------------------------------------- ; May 98: Procedures for reading the SSJ4/SSJ5, SSIES and SSM data were orginally ; coded as FORTRAN subroutines which were linked to the IDL code with ; the call_external function. The FORTRAN code has been replace with ; IDL procedures. ;----------------------------------------------------------------------- ; ; ; PRO particle - procedure to read SSJ4 or SSJ5/Mode A data ; ; Inputs: ; j4file - (string) path/name of SSJ4 or SSJ5 input file ; begtime - start time = Year, Month, Day, Hour, Minute, Julian Day ; totlen = float(endtime - currtime) (minutes) ; x_length - minutes of data to be plotted on each page. May be greater ; than begtime + totlen ; imulf - this is a multiplicative factor, either 1, 2 or 3. ; When the user enters a minlat, sometimes the minutes read ; do not come out as a factor of 2( for time periods be- ; tween 11 and 20) or 3 (for time periods between 21 and ; 30). After you "trip" past the minlat, this factor is ; checked. Say it comes out to 25 minutes read. The logic ; will force the code to read 2 more minutes until 27 ; (imulf*3) are read. This whole set-up was devised to ; force the tick marks to line up with the ehemeris. ; satel - DMSP satellite number ( 6 to 16 or most recent spacecraft) ; mlat - Corr. Mag. Latitude, after tracing along field line to 110 km altitude, ; at start of each minute (limit of 31 minutes) ; mlt - Corr. Mag. Local Time at start of each minute (limit of 31 minutes) ; j4op - flag for opening SSJ4 or SSJ5 file: if j4op = 0, then the file will be ; opened and j4op will be set to 1 ; j4cl - flag for closing SSJ4 or SSJ5 file: if j4cl = 1, then the file will be ; closed and no data retrieved ; minimum_cglat - minimum mag.lat. to be used for plotting ; nsq - flag to indicate which hemisphere to get data from: ; 1=Northern Hemisphere, -1=Southern Hemisphere, 0=Both ; rcd - keeps count of the j4 records which are read in ; iyear - Year of the data file (used to determine geometric factor degraded by ; time on orbit) ; ; OUTPUTS: ; eflux - diff number flux from each of the electron channels ; ecnt - count rate from each of the electron channels ; iflux - diff number flux from each of the ion channels ; icnt - cout rate from each of the ion channels ; plotparsw if 0, no data was found in the minlat range requested ; coinciding with the number of minutes requested ; if 1, the logic stepped thru accumulating eflux, iflux ; and their counts, and icount is > 0 ; removebackground - flag to indicate whether background should be removed ; (NE 0) or not (=0) ;----------------------------------------------------------------------------- pro particle, j4file$, begtime, x_length, imulf, eflux, ecnt, iflux, $ icnt, satel, mlat, mlt, j4op, j4cl, $ ; 10/08 Emery: add more ephemeris data ; minimum_cglat, nsq, totlen, plotparsw, removebackground, iyear minimum_cglat, nsq, totlen, plotparsw, removebackground, iyear, $ dglat,dglon,salt,dglat110,dglon110,dmlon110,uts ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * bdata = intarr(2640) raw2count = fltarr(512) counts = fltarr(40,60) e_energy = fltarr(20) e_geomf = fltarr(20) e_intenergy = fltarr(20) e_intnumber = fltarr(20) e_difnumber = fltarr(20) i_energy = fltarr(20) i_geomf = fltarr(20) i_intenergy = fltarr(20) i_intnumber = fltarr(20) i_difnumber = fltarr(20) jcntin = fltarr(20) jcntout = fltarr(20) abkg = fltarr(2) mode_flag = fltarr(60) modeb_cycle = fltarr(60) num_min_saved_for_plot = 0 modeb_zone = 0 ; - Number of last J5 Mode B sector read out. Set to zero for J4 ; or J5 in Mode A. ;open the J4/J5 data file for input GET_LUN, U openr,U,j4file$,ERROR=err if err ne 0 then begin print,"Cannot open SSJ4/SSJ5 file ",j4file$ goto, done_j4 endif adata =assoc(U, intarr(2640)) lastrecread = 0 ; Create the conversion factors to convert from raw data to counts create_raw2count,raw2count ; set up 'null' values in variables for holding info needed to ; skip records to get to the record wanted. hour = -1 minute = -1 istartminofday = 60*fix(begtime[3]) + fix(begtime[4]) iendminofday = istartminofday + fix(totlen) isatnumber = fix(satel) if isatnumber le 15 then stj4j5 = 'saving J4' else stj4j5 = 'saving J5' ; Get the calibration factors to convert from counts to fluxes read_geomfac, isatnumber, e_energy, e_intenergy, e_intnumber, e_difnumber, $ e_geomf, i_energy, i_intenergy, i_intnumber, i_difnumber, $ i_geomf,iyear ; Set up couters icycle = 0 icycle_delta = 0 idelta_rec = 9999 iskiprecord = 0 istartminofplot = -1 abkg[0:1] = 0 ; Cancel values for previous background computations ; 10/08: TEMP write print,FORMAT='(a10,i5,i3,i7,90f12.3)',"particle: ",iyear,nsq,istartminofday,satel,minimum_cglat, $ x_length,totlen,begtime ; begin processing data. This is the top of a loop. READ_J4_AGAIN: ; read one record of data if not eof(U) then begin on_ioerror,IOERROR_J4READ bdata = adata(icycle) BYTEORDER,BDATA,/NTOHS ;SWAP BYTES 2X2 endif else begin print, "J4/J5: End of File encountered" goto,done_j4 endelse ; convert ephemeris data raw2ephem,bdata,jday,hour,minute,isecond,iyear,glat,glon,alt,gglat110,$ ; 10/08 Emery: add cglon110 ; gglon110,cglat110,mlthr,mltmin,mltsec gglon110,cglat110,cglon110,mlthr,mltmin,mltsec minofday = 60*hour + minute ; 10/08 Emery: save ephemeris data no matter what if (minofday ge istartminofday) and (minofday le istartminofday+x_length) then begin num_min = (minofday - istartminofday) if (num_min gt 30) then print,format='(10i5)',num_min,minofday,istartminofday,iendminofday,begtime mlat[num_min] = cglat110 mlt[num_min] = mlthr+(mltmin+(mltsec/60.))/60. dglat[num_min] = glat dglon[num_min] = glon salt[num_min] = alt dglat110[num_min] = gglat110 dglon110[num_min] = gglon110 dmlon110[num_min] = cglon110 uts[num_min] = float(hour)*3600. + float(minute)*60. + float(isecond) ;print,format='(i5,f8.0)',num_min,uts[num_min] endif ; 10/08: TEMP PRINT ;print,FORMAT='(a6,2i6,2i3,2i5,7f12.3,3i3,f12.3)',"Read: ",icycle,idelta_rec,hour,minute,isecond,iyear, $ ; glat,glon,alt,gglat110,gglon110,cglat110,cglon110,mlthr,mltmin,mltsec,mlthr+(mltmin+(mltsec/60.))/60. ;print,FORMAT='(a6,2i6,2i3)',"Read: ",icycle,idelta_rec,hour,minute if icycle eq 0 then begin ; check the date of the file if jday ne begtime(5) or iyear ne begtime(0) then begin print,'Request date and file date disagree:' print,'Request date =',begtime print,'File date = ', iyear,hour,minute,isecond,jday goto, done_j4 endif endif ; 10/08: TEMP print ;print,FORMAT='(7i6)',minofday,istartminofday,istartminofplot,cglat110,iendminofday, $ ; num_min_saved_for_plot,idelta_rec ; 10/08 Emery: moved after check of idelta_rec=9999 ; Check for start of plot interval ;if (minofday lt istartminofday) or (abs(cglat110) lt minimum_cglat) or $ ; (nsq*cglat110 lt 0) then goto, SKIPMINUTE ; 10/08 Emery: changed gt to ge for idelta_rec=9999 ;if minofday gt istartminofday then begin if minofday ge istartminofday then begin if idelta_rec gt 1 and idelta_rec lt 9999 then begin ;There might be a data gap. If the program has been skipping records. ;Re-set record position to last read before this point and read ;forward without skipping idelta_rec = fix(idelta_rec/2) if idelta_rec lt 1 then idelta_rec = 1 icycle = lastrecread+idelta_rec goto,READ_J4_AGAIN endif idelta_rec = 1 endif ; 10/08: TEMP print ;print,FORMAT='(7i6)',minofday,istartminofday,istartminofplot,cglat110,iendminofday, $ ; num_min_saved_for_plot,idelta_rec ; 10/08 Emery: moved after check of idelta_rec=9999, and allow 3 more mlat of data to get partial mins ; Check for start of plot interval if (minofday lt istartminofday) or (abs(cglat110) lt minimum_cglat-3) or $ (nsq*cglat110 lt 0) then goto, SKIPMINUTE ; Check for end of plot interval if minofday gt iendminofday then goto,done_j4 ; Check for missing minutes of data if num_min_saved_for_plot ne (minofday - istartminofday) then num_min_saved_for_plot = (minofday - istartminofday) if istartminofplot ge 0 then begin ; If next minute of data would be outside the plot frame, then plot the data ; already saved if minofday ge istartminofplot+x_length then goto,done_j4 endif else begin ; If the first minute of data whould be outside the plot frame, then there is nothing to plot if num_min_saved_for_plot gt x_length then begin num_min_saved_for_plot = 0 idelta_rec = 1 goto,done_j4 endif endelse ; reset start time to start of actual data if num_min_saved_for_plot eq 0 then begin begtime[3] = hour begtime[4] = minute endif ; Save ephemeris values mlat[num_min_saved_for_plot] = cglat110 mlt[num_min_saved_for_plot] = mlthr+(mltmin+(mltsec/60.))/60. ; 07/10 skip print ;print,format='(a9,2i4,2f10.1)', stj4j5, hour, minute, mlat(num_min_saved_for_plot),mlt(num_min_saved_for_plot) if num_min_saved_for_plot eq x_length then goto,done_j4 ; Convert minute of data from raw values to counts counts[0:39,0:59] = 0. convert2counts, bdata, raw2count, counts, isatnumber, mode_flag, modeb_cycle,$ modeb_zone ; ; Remove background counts from the J4/J5 data ; if removebackground NE 0 then begin for n=0,59 do begin for itype=0,1 do begin ; Transfer Electron or Ion Counts to temporary array if itype eq 0 then jcntin[0:19] = counts[0:19,n] else $ jcntin[0:19] = counts[20:39,n] ; Compute background and remove it from spectrum. jbackgroundremoval,jcntin,jcntout,isatnumber,abkg,itype ; Transfer Electron or Ion Counts after background removal from temporary array if itype eq 0 then counts[0:19,n] = jcntout[0:19] else $ counts[20:39,n] = jcntout[0:19] endfor ;print,hour,minute,n,abkg endfor endif ; ; store the counts for i=0,59 do begin for j=0,19 do begin ecnt(60*num_min_saved_for_plot+i,j) = counts(j,i) icnt(60*num_min_saved_for_plot+i,j) = counts(20+j,i) endfor endfor ; Convert counts to flux and Save minute of data for future plotting convert2flux, counts, num_min_saved_for_plot, $ eflux, iflux, e_geomf, i_geomf, isatnumber if num_min_saved_for_plot eq 0 then istartminofplot = minofday num_min_saved_for_plot =num_min_saved_for_plot + 1 ; Go read the next minute of data lastrecread = icycle ; 10/08 Emery: change lt to le ;if minofday lt iendminofday-1 then begin if minofday le iendminofday-1 then begin idelta_rec = 1 icycle = icycle + idelta_rec ; increment index to read the next minute of data from file goto,READ_J4_AGAIN endif else goto,done_j4 SKIPMINUTE: ; we got here either becuase data is a low latitude or record ; was earlier than the start of the range if (60*hour+minute) lt istartminofday then begin ; have not reached the start of the data ; Skip ahead and ; read the requested starting time. lastrecread = icycle if idelta_rec gt 1 then begin if isatnumber ge 8 then $ i = fix((istartminofday - minofday)/2) else i = 1 if i lt idelta_rec then idelta_rec = i else idelta_rec = idelta_rec - 1 endif else idelta_rec = 1 if idelta_rec lt 1 then idelta_rec = 1 icycle = icycle + idelta_rec goto,READ_J4_AGAIN endif else begin if minofday gt iendminofday then goto,done_j4 if isatnumber ge 8 and idelta_rec gt 1 then begin ; skip over some of the mid-latitude data icycle_delta = fix((minimum_cglat-abs(cglat110))/3.5) if icycle_delta lt 1 then icycle_delta = 1 endif else begin ; For F6 and F7 data, the mid-latitude data is missing from the file. icycle_delta = 1 endelse if icycle_delta lt 1 then icycle_delta=1 icycle = icycle + idelta_rec goto,READ_J4_AGAIN endelse goto,done_j4 IOERROR_J4READ: ; Becuase of logic to skip records, it is possible to try to read a record that is beyond ; the end-of-file. This causes and IOerror condition. Since error may be solved by ; skipping fewer records, the program needs to clear the IOerror condition flag. This is ; done by closing and re-opening the IO unit. close,U openr,U,j4file$ adata =assoc(U, intarr(2640)) done_j4: if idelta_rec gt 1 and idelta_rec lt 9999 then begin ; If the logic got to this point while trying to skip forward in reading records, then ; cut the size of the number of records skipped and start reading again from the last ; known good record. idelta_rec = fix(idelta_rec/2) if idelta_rec lt 1 then idelta_rec = 1 icycle = lastrecread+idelta_rec goto,READ_J4_AGAIN endif if num_min_saved_for_plot le 0 then plotparsw = 0 else plotparsw = 1 close,U FREE_LUN,U end ;------------------------------------------------------------------------------ ; PRO create_raw2count - create table of conversion factors from SSJ4 raw data ; to counts. Remember that reported count rates greater ; 1.0E6 may be lower than actual counts due to aliasing. ;------------------------------------------------------------------------------ PRO create_raw2count,raw2count for i=0,511 do begin iy = long(i/32) ix = float(i) - iy*32. raw2count(i) = float(ix+32.)*float(2.^iy)-33. endfor end ;------------------------------------------------------------------------------ ; PRO read_geomfac - read the geometric factors for the SSJ4/SSJ5 instrument. ; These factors are used to convert counts to differential ; fluxes. For SSJ5, the factors are only for Mode A. These ; factors used to be in a auxiliary file; now they are ; coded into a sub-procedure (getjgf). ;------------------------------------------------------------------------------ PRO read_geomfac, isatnumber, e_energy, e_intenergy, e_intnumber, e_difnumber, $ e_geomf, i_energy, i_intenergy, i_intnumber, i_difnumber, $ i_geomf, iyear, JGF_LABEL ;------------------------------------------------------------------------------ ; isatnumber - DMSP spacecraft number. (for DMSP F10, isatnumber = 10) ; e_energy - Energy of the 20 Electron Channels (keV) ; i_energy - Energy of the 20 Ion Channels (keV) ; e_geomf - geometric factors for the SSJ4 and SSJ5/Mode A electron data (20 channels) ; i_geomf_b - geometric factors for the SSJ4 and SSJ5/Mode A ion data (20 channels) ; e_intenergy - conversion factor from counts to integral electron energy flux ; e_intnumber - conversion factor from counts to integral electron number flux ; e_difnumber - conversion factor to differential electron number flux ; i_intenergy - conversion factor to integral ion number flux ; i_intnumber - conversion factor from counts to integral ion number flux ; i_difnumber - conversion factor to differential ion number flux ; iyear - year of the data ; JGF_LABEL - text string indicating the date that the SSJ calibration was released ; ; Comments: ; Feb 2004 - Changed this routine to get values from procedure getjgf instead ; of reading a file called 'geomfac.dat ;------------------------------------------------------------------------------ GME = fltarr(22) GMI = fltarr(22) EE = fltarr(22) EI = fltarr(22) ZF = fltarr(6,2) engw = fltarr(20) engi = fltarr(20) CF = fltarr(4) iflg = 0 iflg = 129 GETJGF,ISATnumber,iyear,GME,GMI,EE,EI,IFLG,ZF,CF,JGF_LABEL for i=0,9 do begin e_energy[i] = ee[i+3] i_energy[i] = ei[i+3] e_geomf[i] = gme[i+3]*1.e-3 i_geomf[i] = gmi[i+3]*1.e-3 endfor e_energy[10] = e_energy[9] i_energy[10] = i_energy[9] e_geomf[10] = 0. i_geomf[10] = 0. for i=11,19 do begin e_energy[i] = ee[i+2] i_energy[i] = ei[i+2] e_geomf[i] = gme[i+2]*1.e-3 i_geomf[i] = gmi[i+2]*1.e-3 endfor if isatnumber lt 16 then sampletime = 0.098 else sampletime = 0.050 e_CONST=1.603E-37 i_CONST=5.399E-31 for i=1,18 do begin engw[i] = (e_energy[i-1]-e_energy[i+1])/2. engi[i] = (i_energy[i-1]-i_energy[i+1])/2 endfor engw[0] = e_energy[0]-e_energy[1] engw[19] = e_energy[18]-e_energy[19] engw[9] = (e_energy[8]-e_energy[11])/2. engw[10] = engw[9] engi[0] = i_energy[0]-i_energy[1] engi[19] = i_energy[18]-i_energy[19] engi[9] = (i_energy[8]-i_energy[11])/2. engi[10] = engw[9] for i=0,19 do begin if e_geomf[i] gt 0. then e_difnumber[i] = 1./(e_geomf[i]*sampletime) if i_geomf[i] gt 0. then i_difnumber[i] = 1./(i_geomf[i]*sampletime) if e_geomf[i] gt 0. then e_intnumber[i] = e_difnumber[i]*engw[i] if e_geomf[i] gt 0. then e_intenergy[i] = e_intnumber[i]*e_energy[i] if i_geomf[i] gt 0. then i_intnumber[i] = i_difnumber[i]*engi[i] if i_geomf[i] gt 0. then i_intenergy[i] = i_intnumber[i]*i_energy[i] endfor end PRO GETJGF,IDSAT,NYR,GME,GMI,EE,EI,IFLG,ZF,CF,JGF_LABEL ;*********************************************************************** ; NOTES: ; ; CHANNEL ENERGIES AND GEOMETRIC FACTORS ARE SHIFTED RIGHT THREE ; CHANNELS IN ANTICIPATION OF OPTIONAL HIGH ENERGY EXTRAPOLATION ; WHEN HIGH ENERGY FLUXES ARE SIGNIFICANT. CHANNEL ELEVEN IS ; DISCARDED AS THE SIMPLEST SOLUTION TO THE J4 OVERLAP CHANNEL AND ; TO BE CONSISTENT WITH THE NINETEEN CHANNEL J5. ; ; THE CORY ARRAY IS THE FINAL RESULTS OF A STUDY OF THE LONG TERM ; VARIATIONS OF THE COUNT RATES IN THE SSJ4 SENSORS. THESE LONG TERM ; VARIATIONS ARE DUE TO VARIATIONS IN THE MEV PARTICLES IN THE INNER ; RADIATION BELTS AND DUE TO DETERIATION WITH TIME OF THE SENSITIVITY ; OF THE FOUR DETECTORS OF EACH J4 INSTRUMENT. THE YEARLY AVERAGE ; OF A COUNT RATES IN THE SOUTH ATLANTIC ANOMALY ARE NORMALIZED TO ; THE HIGH ENERGY ION HEAD AND THEN TO THE FIRST YEAR OF OPERATION ; FOR EACH INSTRUMENT. THE CORRECTION APPLIED TO THE GEOMETRIC ; FACTORS ASSUMES THAT THE PUBLISHED GEOMETRIC FACTORS APPLY TO THE ; FIRST YEAR'S DATA AND THAT A CORRECTION IS JUSTIFIED WHEN THE YEARLY ; INDEX FALLS BELOW ONE. ; ; THIS FILE CONTAINS A REVISION OF THE GEOMETRIC FACTORS FOR ALL THE ; SSJ4 INSTRUMENTS. THIS REVISION WAS DUE TO A CAREFUL RE-EXAMINATION ; ON THE CALIBRATION DATA. THERE WERE SOME SMALL TO MODERATE ERRORS ; MADE IN THE ORIGINAL CREATION OF THE GEOMETRIC FACTOR FROM TEST ; CHAMBER DATA. THE BIGGEST CHANGES ARE IN THE LOWEST ENERGY CHANNELS ; ; FEBRUARY 2004 ;--------- ; ; UPDATE 29 NOV 05 UPDATED CORY FOR 2004 AND 2005 AND UPDATED F16 DATA ; FORTRAN VERSION ;***** ; UPDATE 15 DEC 05 ADDED DATA FOR F17 AND ADJUSTED LOGIC ACCORDINGLY ; FORTRAN VERSION F17 DATA WILL HAVE TO BE ADJUSTED AFTER LAUNCH ;***** ; UPDATE 19 DEC 05 ADDED YEAR 2006 TO ESTIMATE TO CORY BASED ON 2004, ; FORTRAN VERSION 2005 BEHAVIOR ; ; 13 APRIL 2006 UPDATED IDL VERSION ;-------- ; UPDATE 15 MAY 06 UPDATED F16, F17, F20 GFs AND ZPs ; ------- ; UPDATE NOV 06 ADDED (GFL) INFLIGHT CALIBRATION VALUES BASED ON ; HARDY PROJECT ANALYSIS. THEY ARE APPLIED ONLY TO THE ; LOW ENERGY CHANNELS BUT ARE PROBABLY APPROPRIATE FOR ; ALL ; ------- ; UPDATE 14 MAR 07 FINAL 2006 CORY VALUES, FINAL CALIBRATION FOR F17 ; AND PRELIMINARY CALIBRATION FOR F18 ADDED ; ------- ; UPDATE 08 JAN 2008 UPDATED BASED ON 14 DEC 07 FORTRAN UPDATE. CORY ; VALUES BASED ON FURTHER ANALYSIS USING NEW HARDY ; MODEL DATABASE AND THE NEW UNDERSTANDING THAT THE ; SSJ5 MICROCHANNEL PLATES DEGRADE IN A ENERGY ; DEPENDENT MANNER ; F12 CORY VALUES REPLACED WITH SET DERIVED FROM SOLAR ; CYCLE VARIATION SEEN IN ANOMALY INDEX. ; ADDED HIGH ENERGY ION CORY VALUES DERIVED FROM SSN ; VARIATION. ; 'RAMP' TYPE CORY DEFINED FOR J5 IONS AND ELECTRONS. ; HIGH ENERGY ION AND ELECTRON CORY REPLACED BY ; 0.75 + CORY / 4 (AVERAGE RAMP VALUE). ; IFC (XEIFC) REDERIVED AND APPLIED TO ALL CHANNELS ; PRELIMINARY IFC (XIIFC) ADDED FROM IONS ; ------- ; UPDATE 05 MAY 08 FINAL 2007 CORY VALUES and PRELIMINARY 2008 ; ADDED F18 WITH PRELIMINARY GFs ; ; UPDATE 22 SEP 08 ADDED PRELIMINARY 2008 (8 MONTH) CORY ; ; UPDATE 03 FEB 09 FINAL 2008 CORY VALUES ; ; FREDERICK RICH, AFRL, HANSCOM AFB, MA (frederick.rich@hanscom.af.mil) ; ERNEST HOLMAN, BOSTON COLLEGE, CHESTNUT HILL, MA ; (ernest.holman.ctr@hanscom.af.mil) ; KATHARINE KADINSKY-CADE, AFRL, HANSCOM AFB, MA ; (katharine.kadinsky-cade@hanscom.af.mil) ;***** ; INPUT ARGUMENTS ; IDSAT SATELLITE ID: MAY BE 'F' NUMBER IN RANGE 6 <= IDSAT <= 18 ; OR 4-DIGIT DESIGNATOR ; NYR YEAR NUMBER: EITHER 2-DIGIT OR 4-DIGIT ; IFLG CONTROL FLAG BIT1 SET IMPLIES DO CORY CORRECTIONs ; BIT8 SET IMPLIES DO LOW ENERGY GF CORRECTION ; ; OUTPUT ARRAYS ; GME ELECTRON CHANNEL GEOMETRIC FACTORS ; GMI ION CHANNEL GEOMETRIC FACTORS ; EE ELECTRON CHANNEL ENERGIES ; EI ION CHANNEL ENERGIES ; ZF ZONE SCALE FACTORS FOR J5 INSTRUMENTS ; CF CORRECTION FACTORS APPLIED TO EACH OF FOUR HEADS ; (not used in March 2007 update in order to maintain ; compatibility with older versions but is used in ; Dec 07-Jan08 update) JGF_LABEL = 'Cal as of Dec 07' ;*********************************************************************** ; ZF = fltarr(6,2) ; GME = fltarr(22) ; GMI = fltarr(22) ; EE = fltarr(22) ; EI = fltarr(22) ; CF = fltarr(4) GG = fltarr(20,2,13) CORY = fltarr(4,23,13) ZFAC = fltarr(6,2,3) ; Energy (keV) of Channel No. First three channels don't exist in the SSJ4/J5 data ; but are included for the purpose of extrapolating the data to higher energy. EGY = [95.0,65.0,44.0,30.0,20.4,13.9,9.45,6.46,4.40,3.00,2.04,1.392,.949,.646,.440,.300,.204,.139,.095,.065,.044,.030] ;;GFL = [1.418,.846,1.128,1.354,.940,.540,.623,1.131,1.354,1.542,1.0,1.0,1.0] ; ; set of in-flight-calibration parameters for electrons derived from the Hardy Project average spectra XEIFC = [1.409,1.182,0.845,1.230,0.632,0.522,0.743,1.000,1.302,1.237,1.092,1.191,1.000] ; set of in-flight-calibration parameters for ions derived from the Hardy Project average spectra XIIFC = [0.724,0.658,0.552,.516,0.681,0.881,0.850,1.642,0.767,2.081,1.000,1.286,1.000] ; ; 4-DIGIT SATELLITE DESIGNATORS IDS = [7540,8541,9543,0542,1544,2546,3545,4547,5548,6549,7554,8551,9551] ; F18 designators is a guess ; ; F06 GG[0:19,0,0]= [0.75,0.49,0.41,0.33,0.27,0.21,0.16,0.13,0.096,0.076,0.0157,.0118,.009423,.006121,.004153,.002019,.0007967,.0003058,.00008534,.00001523] GG[0:19,1,0]= [1.8,1.25,.8439,.5730,.3959,.2709,.1875,.1250,.08439,.05835,2.437,1.589,1.483,.8369,.5827,.3920,.2649,.1907,.1271,.08475] ; ; F07 GG[0:19,0,1] = [ 0.58,0.49,0.41,0.33,0.27,0.21,0.16,0.13,0.096,0.076,.032,.0255,.02036,.01349,.009610,.005061,.002342,.001058,.0004058,.0001191] GG[0:19,1,1] = [ 2.4,1.667,1.146,.7814,.5418,.3647,.2501,.1771,.1146,.07918,2.058,1.372,.9603,.6467,.4507,.3038,.2058,.1470,.09799,.06859] ; ; F08 GG[0:19,0,2] = [ 0.326,0.275,0.23,0.185,0.152,0.118,0.0899,0.073,0.0539,0.0427,.03163,.02526,.02013,.01332,.009514,.005007,.002316,.001046,.0004023,.0001178] GG[0:19,1,2] = [ 1.15,0.8012,.5512,.3761,.2605,.1750,.1198,.08512,.05512,.03803,1.096,.7315,.5123,.3445,.2399,.1619,.1096,.07838,.05222,.03653] ; ; F09 GG[0:19,0,3] = [ 0.201,0.17,0.142,0.115,0.0937,0.0729,0.0555,0.0451,.0333,.0264,.0369,.0294,.0235,.0156,.0110,.00583,.00270,.00122,.000468,.000137] GG[0:19,1,3] = [ 1.14,0.7939,.5459,.3720,.2584,.1740,.1188,.08439,.05459,.03772,1.266,.8411,.5890,.3962,.2765,.1859,.1266,.09015,.06007,.04206] ; ; F10 GG[0:19,0,4] = [ .462,.389,.302,.243,.194,.158,.126,.102,.0824,.0669,.05206,.03714,.02506,.01740,.01219,.006839,.003697,.001952,.0009389,.0003628] GG[0:19,1,4] = [ .5670,.4032,.2709,.1823,.1240,.08523,.05772,.03938,.02678,.01844,.5549,.3755,.2548,.1739,.1196,.08097,.05549,.03821,.02636,.01794] ; ; F11 GG[0:19,0,5] = [ .631,.544,.428,.349,.278,.226,.182,.149,.121,.0950,.1092,.07788,.05345,.03585,.02460,.01391,.007392,.003907,.001846,.0006987] GG[0:19,1,5] = [ .5490,.3970,.2709,.1844,.1261,.08585,.05876,.04032,.02771,.01844,.5818,.4065,.2787,.1939,.1324,.09093,.06120,.04308,.02973,.02056] ; ; F12 GG[0:19,0,6] = [ .4552,.3815,.3188,.2661,.2232,.1858,.1561,.1297,.1088,.09103,.06909,.04888,.03279,.02248,.01536,.008516,.004523,.002338,.001094,.0004161] GG[0:19,1,6] = [ .7060,.5022,.3428,.2344,.1594,.1094,.07449,.05084,.03469,.02375,.6390,.4371,.2988,.2051,.1402,.09600,.06566,.04499,.03084,.02098] ; ; F13 GG[0:19,0,7] = [ .456,.362,.287,.228,.181,.144,.114,.0908,.0721,.0572,.03644,.02615,.01777,.01232,.008563,.004827,.002598,.001359,.0006463,.0002496] GG[0:19,1,7] = [ .984,.6981,.4761,.3240,.2209,.1511,.1027,.07001,.04772,.03251,.05572,.03798,.02602,.01778,.01214,.008305,.005679,.003869,.002651,.001810] ; ; F14 GG[0:19,0,8] = [ 0.334,0.272,0.222,0.181,0.148,0.121,0.0985,0.0803,0.0656,0.0535,.03939,.02828,.01930,.01345,.009314,.005287,.002839,.001488,.0007104,.0002746] GG[0:19,1,8] = [ 1.33,.9408,.6397,.4355,.2959,.2011,.1365,.09314,.06335,.04303,1.008,.6877,.4705,.3206,.2196,.1499,.1021,.06970,.04763,.03253] ; ; F15 GG[0:19,0,9] = [ 0.3124,0.2528,0.2046,0.1657,0.1341,0.1086,0.08790,0.07116,0.0576,0.04662,.05072,.03647,.02488,.01733,.01205,.006807,.003672,.001928,.0009189,.0003560] GG[0:19,1,9] = [ 0.9016,0.6392,.435,.296,.2014,.1371,.09328,.06347,.04320,.02939,0.9439,0.6443,0.4398,0.3002,0.2049,0.1398,0.09454,0.06515,0.04447,0.03035] ; F16 (J5S16) ;GG[0:19,0,10] = [1.400,1.161,0.934,0.735,0.567,0.433,0.327,0.241,0.177,0.131,9.999,0.0968,0.0688,0.0482,0.0337,0.0227,0.0143,0.00889,0.00488,0.00241] ; Prior to March 2007 ;GG[0:19,1,10] = [ 18.9,12.1,7.66,4.80,3.00,1.89,1.20,0.768,0.503,0.349, 9.999, 0.248, 0.166, 0.117, 0.0800, 0.0482, 0.0278, 0.0168, 0.00787, 0.00458] GG[0:19,0,10] = [1.993,1.653,1.330,1.047,0.808,0.617,0.465,0.343,0.252,0.186,9.999,0.138,0.0980,0.0686,0.0480,0.0323,0.0204,0.0127,0.00695,0.00343] GG[0:19,1,10] = [ 21.7, 13.9, 8.86, 5.60, 3.58, 2.29, 1.47,0.938,0.601,0.398,9.999,0.264, 0.171, 0.117,0.0824,0.0559,0.0376,0.0257, 0.0122,0.00643] ; F17 (J5S20) ;GG[0:19,0,11] = [0.540,0.407,0.307,0.234,0.196,0.148,0.111,0.0782,0.0588,0.0426,9.999,0.0305,0.0218,0.0149,0.0102,0.00752,0.00517,0.00403,0.00273,0.00156] ;GG[0:19,1,11] = [10.6,6.97,4.57,2.99,1.92,1.28,0.832,0.525,0.326,0.224,9.999,0.155,0.104,0.0713,0.0457,0.0292,0.0178,0.00745,0.00277,0.00105] GG[0:19,0,11] = [1.177,0.911,0.706,0.539,0.411,0.308,0.225,0.168,0.121,0.0847,9.999,0.0635,0.0463,0.0324,0.0220,0.0157, 0.0106,0.00760,0.00528, 0.00349] GG[0:19,1,11] = [ 9.50, 6.33, 4.22, 2.82, 1.87, 1.18,0.777,0.507,0.330, 0.201,9.999, 0.149,0.0962,0.0511,0.0285,0.0161,0.00843,0.00515,0.00272,0.000996] ; F18 (J5S17) GG[0:19,0,12] = [0.997,0.742,0.570,0.433,0.367,0.270,0.206,0.149,0.110,0.0789,9.999,0.0566,0.0408,0.0279,0.0192,0.0143,0.00974,0.00773,0.00531,0.00328] GG[0:19,1,12] = [ 10.3, 6.71, 4.38, 2.73, 1.76, 1.15,0.752,0.471,0.288, 0.202,9.999, 0.146, 0.102,0.0705,0.0435,0.0315, 0.0209, 0.0109,0.00435,0.00177] ; ; MODE B ZONE PROFILES FOR J5s ;ZFAC = [0.716,0.828,1.004,1.239,1.330,1.189,0.942,1.112,1.178,1.070,1.009,0.791] ZFAC[0:5,0,0] = [ 0.661,0.734,0.864,1.089,1.306,1.346] ZFAC[0:5,1,0] = [ 1.012,1.134,1.040,0.970,0.991,0.853] ZFAC[0:5,0,1] = [ 0.976,0.790,0.667,0.915,1.184,1.469] ZFAC[0:5,1,1] = [ 0.935,1.042,0.991,1.012,1.003,1.017] ZFAC[0:5,0,2] = [ 1.898,1.198,0.784,0.679,0.664,0.777] ZFAC[0:5,1,2] = [ 1.097,1.094,0.970,0.883,0.971,0.985] ; ; CORG MATRICES FOR F6 THROUGH F18 CORY[0:3, 0, 0] = [1.000,1.000,1.000,1.000] ; F6 1986 ;CORY[0:3, 1, 0] = [1.000,0.995,1.000,0.998] ; 1987 CORY[0:3, 1, 0] = [1.000,0.995,0.981,0.998] ; 1987 CORY[0:3, 2, 0] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3, 0] = [0.000,0.000,0.000,0.000] CORY[0:3, 4, 0] = [0.000,0.000,0.000,0.000] CORY[0:3, 5, 0] = [0.000,0.000,0.000,0.000] CORY[0:3, 6, 0] = [0.000,0.000,0.000,0.000] CORY[0:3, 7, 0] = [0.000,0.000,0.000,0.000] CORY[0:3, 8, 0] = [0.000,0.000,0.000,0.000] CORY[0:3, 9, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,10, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,11, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,12, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,13, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,14, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,15, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,16, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,17, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,18, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,19, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,20, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,21, 0] = [0.000,0.000,0.000,0.000] CORY[0:3,22, 0] = [0.000,0.000,0.000,0.000] CORY[0:3, 0, 1] = [1.000,1.000,1.000,1.000] ; F7 1986 ;CORY[0:3, 1, 1] = [1.000,1.000,1.000,1.000] ; 1987 ;CORY[0:3, 2, 1] = [0.961,1.000,1.000,0.985] ; 1988 CORY[0:3, 1, 1] = [1.000,1.000,0.973,1.000] ; 1987 CORY[0:3, 2, 1] = [0.988,1.000,1.000,0.985] ; 1988 CORY[0:3, 3, 1] = [0.000,0.000,0.000,0.000] CORY[0:3, 4, 1] = [0.000,0.000,0.000,0.000] CORY[0:3, 5, 1] = [0.000,0.000,0.000,0.000] CORY[0:3, 6, 1] = [0.000,0.000,0.000,0.000] CORY[0:3, 7, 1] = [0.000,0.000,0.000,0.000] CORY[0:3, 8, 1] = [0.000,0.000,0.000,0.000] CORY[0:3, 9, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,10, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,11, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,12, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,13, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,14, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,15, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,16, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,17, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,18, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,19, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,20, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,21, 1] = [0.000,0.000,0.000,0.000] CORY[0:3,22, 1] = [0.000,0.000,0.000,0.000] CORY[0:3, 0, 2] = [0.000,0.000,0.000,0.000] ; F8 CORY[0:3, 1, 2] = [1.000,1.000,1.000,1.000] ; 1987 ;CORY[0:3, 2, 2] = [0.884,1.000,1.000,1.000] ; 1988 ;CORY[0:3, 3, 2] = [0.970,1.000,1.000,1.000] ; 1989 ;CORY[0:3, 4, 2] = [0.997,1.000,1.000,1.000] ; 1990 ;CORY[0:3, 5, 2] = [0.999,0.973,1.000,1.000] ; 1991 ;CORY[0:3, 6, 2] = [0.943,0.858,1.000,1.000] ; 1992 ;CORY[0:3, 7, 2] = [0.937,0.857,1.000,1.000] ; 1993 ;CORY[0:3, 8, 2] = [0.937,0.864,1.000,1.000] ; 1994 CORY[0:3, 2, 2] = [1.000,1.000,1.000,1.000] ; 1988 CORY[0:3, 3, 2] = [1.000,1.000,0.920,1.000] ; 1989 CORY[0:3, 4, 2] = [1.000,1.000,0.946,1.000] ; 1990 CORY[0:3, 5, 2] = [1.000,1.000,0.946,1.000] ; 1991 CORY[0:3, 6, 2] = [1.000,0.987,0.882,1.000] ; 1992 CORY[0:3, 7, 2] = [1.000,0.987,0.944,1.000] ; 1993 CORY[0:3, 8, 2] = [1.000,0.994,0.891,1.000] ; 1994 CORY[0:3, 9, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,10, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,11, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,12, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,13, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,14, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,15, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,16, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,17, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,18, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,19, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,20, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,21, 2] = [0.000,0.000,0.000,0.000] CORY[0:3,22, 2] = [0.000,0.000,0.000,0.000] CORY[0:3, 0, 3] = [0.000,0.000,0.000,0.000] ; F9 CORY[0:3, 1, 3] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2, 3] = [1.000,1.000,1.000,1.000] ; 1988 ;CORY[0:3, 3, 3] = [1.000,1.000,1.000,1.000] ; 1989 ;CORY[0:3, 4, 3] = [1.000,1.000,1.000,1.000] ; 1990 ;CORY[0:3, 5, 3] = [1.000,0.937,1.000,1.000] ; 1991 ;CORY[0:3, 6, 3] = [1.000,0.927,1.000,1.000] ; 1992 ;CORY[0:3, 7, 3] = [0.937,0.857,1.000,1.000] ; 1993 ;CORY[0:3, 8, 3] = [0.937,0.864,1.000,1.000] ; 1994 CORY[0:3, 3, 3] = [1.000,1.000,0.836,1.000] ; 1989 CORY[0:3, 4, 3] = [1.000,1.000,0.832,1.000] ; 1990 CORY[0:3, 5, 3] = [1.000,0.937,0.873,1.000] ; 1991 CORY[0:3, 6, 3] = [1.000,0.927,0.792,1.000] ; 1992 CORY[0:3, 7, 3] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8, 3] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,10, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,11, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,12, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,13, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,14, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,15, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,16, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,17, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,18, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,19, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,20, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,21, 3] = [0.000,0.000,0.000,0.000] CORY[0:3,22, 3] = [0.000,0.000,0.000,0.000] CORY[0:3, 0, 4] = [0.000,0.000,0.000,0.000] ; F10 CORY[0:3, 1, 4] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2, 4] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3, 4] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4, 4] = [1.000,1.000,1.000,1.000] ; 1990 CORY[0:3, 5, 4] = [1.000,0.618,1.000,1.000] ; 1991 ;CORY[0:3, 6, 4] = [0.956,0.797,1.000,0.967] ; 1992 ;CORY[0:3, 7, 4] = [0.899,0.906,1.000,0.844] ; 1993 ;CORY[0:3, 8, 4] = [0.852,0.811,1.000,0.469] ; 1994 ;CORY[0:3, 9, 4] = [0.854,0.721,1.000,0.179] ; 1995 ;CORY[0:3,10, 4] = [0.829,0.589,1.000,0.067] ; 1996 ;CORY[0:3,11, 4] = [0.823,0.516,1.000,0.016] ; 1997 CORY[0:3, 6, 4] = [0.956,1.000,0.979,0.967] ; 1992 CORY[0:3, 7, 4] = [0.899,1.000,1.029,0.844] ; 1993 CORY[0:3, 8, 4] = [0.852,1.000,1.029,0.469] ; 1994 CORY[0:3, 9, 4] = [0.854,1.000,0.967,0.179] ; 1995 CORY[0:3,10, 4] = [0.829,0.953,0.968,0.067] ; 1996 CORY[0:3,11, 4] = [0.823,0.836,0.968,0.016] ; 1997 CORY[0:3,12, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,13, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,14, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,15, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,16, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,17, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,18, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,19, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,20, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,21, 4] = [0.000,0.000,0.000,0.000] CORY[0:3,22, 4] = [0.000,0.000,0.000,0.000] CORY[0:3, 0, 5] = [0.000,0.000,0.000,0.000] ; F11 CORY[0:3, 1, 5] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2, 5] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3, 5] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4, 5] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5, 5] = [1.000,1.000,1.000,1.000] ; 1991 CORY[0:3, 6, 5] = [1.000,1.000,1.000,1.000] ; 1992 CORY[0:3, 7, 5] = [1.000,1.000,1.000,1.000] ; 1993 ;CORY[0:3, 8, 5] = [0.985,1.000,1.000,1.000] ; 1994 ;CORY[0:3, 9, 5] = [0.938,1.000,1.000,1.000] ; 1995 CORY[0:3, 8, 5] = [0.988,1.000,0.920,1.000] ; 1994 CORY[0:3, 9, 5] = [0.941,1.000,0.847,1.000] ; 1995 CORY[0:3,10, 5] = [0.000,0.000,0.000,0.000] ; 1996 ;CORY[0:3,11, 5] = [0.990,1.000,1.000,0.969] ; 1997 ;CORY[0:3,12, 5] = [0.955,1.000,1.000,0.887] ; 1998 ;CORY[0:3,13, 5] = [0.918,1.000,1.000,0.746] ; 1999 ;CORY[0:3,14, 5] = [1.005,0.962,1.000,0.524] ; 2000 CORY[0:3,11, 5] = [0.993,1.000,0.753,0.990] ; 1997 CORY[0:3,12, 5] = [0.957,1.000,0.728,0.906] ; 1998 CORY[0:3,13, 5] = [0.921,1.000,0.789,0.762] ; 1999 CORY[0:3,14, 5] = [1.000,1.000,0.827,0.536] ; 2000 CORY[0:3,15, 5] = [0.000,0.000,0.000,0.000] CORY[0:3,16, 5] = [0.000,0.000,0.000,0.000] CORY[0:3,17, 5] = [0.000,0.000,0.000,0.000] CORY[0:3,18, 5] = [0.000,0.000,0.000,0.000] CORY[0:3,19, 5] = [0.000,0.000,0.000,0.000] CORY[0:3,20, 5] = [0.000,0.000,0.000,0.000] CORY[0:3,21, 5] = [0.000,0.000,0.000,0.000] CORY[0:3,22, 5] = [0.000,0.000,0.000,0.000] CORY[0:3, 0, 6] = [0.000,0.000,0.000,0.000] ; F12 CORY[0:3, 1, 6] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2, 6] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3, 6] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4, 6] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5, 6] = [0.000,0.000,0.000,0.000] ; 1991 CORY[0:3, 6, 6] = [0.000,0.000,0.000,0.000] ; 1992 CORY[0:3, 7, 6] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8, 6] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9, 6] = [1.000,1.000,1.000,1.000] ; 1995 ;CORY[0:3,10, 6] = [1.000,1.000,1.000,1.000] ; 1996 ;CORY[0:3,11, 6] = [1.000,0.946,1.000,1.000] ; 1997 ;CORY[0:3,12, 6] = [1.000,0.805,1.000,1.000] ; 1998 ;CORY[0:3,13, 6] = [0.988,0.577,1.000,1.000] ; 1999 CORY[0:3,10, 6] = [0.965,0.960,0.933,1.000] ; 1996 CORY[0:3,11, 6] = [0.905,0.840,0.863,1.000] ; 1997 CORY[0:3,12, 6] = [0.853,0.664,0.802,0.974] ; 1998 CORY[0:3,13, 6] = [0.760,0.430,0.724,0.935] ; 1999 CORY[0:3,14, 6] = [0.822,0.295,1.000,1.000] ; 2000 ;CORY[0:3,15, 6] = [0.500,0.103,1.000,1.000] ; 2001 ;CORY[0:3,16, 6] = [0.422,0.060,1.000,1.000] ; 2002 CORY[0:3,15, 6] = [0.526,0.105,0.990,1.000] ; 2001 CORY[0:3,16, 6] = [0.127,0.017,0.284,0.462] ; 2002 CORY[0:3,17, 6] = [0.000,0.000,0.000,0.000] CORY[0:3,18, 6] = [0.000,0.000,0.000,0.000] CORY[0:3,19, 6] = [0.000,0.000,0.000,0.000] CORY[0:3,20, 6] = [0.000,0.000,0.000,0.000] CORY[0:3,21, 6] = [0.000,0.000,0.000,0.000] CORY[0:3,22, 6] = [0.000,0.000,0.000,0.000] CORY[0:3, 0, 7] = [0.000,0.000,0.000,0.000] ; F13 CORY[0:3, 1, 7] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2, 7] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3, 7] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4, 7] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5, 7] = [0.000,0.000,0.000,0.000] ; 1991 CORY[0:3, 6, 7] = [0.000,0.000,0.000,0.000] ; 1992 CORY[0:3, 7, 7] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8, 7] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9, 7] = [1.000,1.000,1.000,1.000] ; 1995 ;CORY[0:3,10, 7] = [1.000,1.000,1.000,1.000] ; 1996 ;CORY[0:3,11, 7] = [1.000,0.982,1.000,0.995] ; 1997 ;CORY[0:3,12, 7] = [0.995,0.905,1.000,0.906] ; 1998 ;CORY[0:3,13, 7] = [0.968,0.839,1.000,0.603] ; 1999 ;CORY[0:3,14, 7] = [0.933,0.764,1.000,0.353] ; 2000 ;CORY[0:3,15, 7] = [0.922,0.738,1.000,0.325] ; 2001 ;CORY[0:3,16, 7] = [0.908,0.686,1.000,0.299] ; 2002 ;CORY[0:3,17, 7] = [0.827,0.593,1.000,0.437] ; 2003 ;CORY[0:3,18, 7] = [0.822,0.614,1.000,0.572] ; 2004 ;CORY[0:3,19, 7] = [0.829,0.617,1.000,0.665] ; 2005 ;CORY[0:3,20, 7] = [0.852,0.685,1.000,0.629] ; 2006 CORY[0:3,10, 7] = [1.000,0.991,0.948,1.000] ; 1996 CORY[0:3,11, 7] = [1.000,0.973,0.962,1.000] ; 1997 CORY[0:3,12, 7] = [1.000,0.898,0.938,0.943] ; 1998 CORY[0:3,13, 7] = [0.988,0.832,0.936,0.628] ; 1999 CORY[0:3,14, 7] = [0.952,0.757,0.944,0.367] ; 2000 CORY[0:3,15, 7] = [0.941,0.732,0.940,0.338] ; 2001 CORY[0:3,16, 7] = [0.927,0.680,0.911,0.311] ; 2002 CORY[0:3,17, 7] = [0.844,0.588,0.928,0.455] ; 2003 CORY[0:3,18, 7] = [0.839,0.609,0.996,0.596] ; 2004 CORY[0:3,19, 7] = [0.845,0.612,0.912,0.698] ; 2005 CORY[0:3,20, 7] = [0.870,0.653,0.872,0.787] ; 2006 CORY[0:3,21, 7] = [0.863,0.644,0.872,0.840] ; 2007 CORY[0:3,22, 7] = [0.854,0.633,0.866,0.885] ; 2008 CORY[0:3, 0, 8] = [0.000,0.000,0.000,0.000] ; F14 CORY[0:3, 1, 8] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2, 8] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3, 8] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4, 8] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5, 8] = [0.000,0.000,0.000,0.000] ; 1991 CORY[0:3, 6, 8] = [0.000,0.000,0.000,0.000] ; 1992 CORY[0:3, 7, 8] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8, 8] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9, 8] = [0.000,0.000,0.000,0.000] ; 1995 CORY[0:3,10, 8] = [0.000,0.000,0.000,0.000] ; 1996 CORY[0:3,11, 8] = [1.000,1.000,1.000,1.000] ; 1997 ;CORY[0:3,12, 8] = [0.987,0.964,1.000,0.963] ; 1998 ;CORY[0:3,13, 8] = [0.974,0.922,1.000,0.933] ; 1999 ;CORY[0:3,14, 8] = [0.957,0.787,1.000,0.932] ; 2000 ;CORY[0:3,15, 8] = [0.941,0.653,1.000,0.870] ; 2001 ;CORY[0:3,16, 8] = [0.944,0.414,1.000,0.670] ; 2002 ;CORY[0:3,17, 8] = [0.866,0.139,1.000,0.612] ; 2003 ;CORY[0:3,18, 8] = [0.856,0.054,1.000,0.579] ; 2004 ;CORY[0:3,19, 8] = [0.940,0.042,1.000,0.530] ; 2005 CORY[0:3,12, 8] = [0.987,0.964,0.974,0.963] ; 1998 CORY[0:3,13, 8] = [0.974,0.922,0.985,0.933] ; 1999 CORY[0:3,14, 8] = [0.957,0.787,0.972,0.932] ; 2000 CORY[0:3,15, 8] = [0.941,0.653,0.970,0.870] ; 2001 CORY[0:3,16, 8] = [0.944,0.414,0.926,0.670] ; 2002 CORY[0:3,17, 8] = [0.866,0.139,0.914,0.612] ; 2003 CORY[0:3,18, 8] = [0.856,0.054,0.945,0.579] ; 2004 CORY[0:3,19, 8] = [0.940,0.042,0.837,0.530] ; 2005 CORY[0:3,20, 8] = [0.000,0.000,0.000,0.000] ; 2006 CORY[0:3,21, 8] = [0.000,0.000,0.000,0.000] ; 2007 CORY[0:3,22, 8] = [0.000,0.000,0.000,0.000] ; 2008 CORY[0:3, 0, 9] = [0.000,0.000,0.000,0.000] ; F15 CORY[0:3, 1, 9] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2, 9] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3, 9] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4, 9] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5, 9] = [0.000,0.000,0.000,0.000] ; 1991 CORY[0:3, 6, 9] = [0.000,0.000,0.000,0.000] ; 1992 CORY[0:3, 7, 9] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8, 9] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9, 9] = [0.000,0.000,0.000,0.000] ; 1995 CORY[0:3,10, 9] = [0.000,0.000,0.000,0.000] ; 1996 CORY[0:3,11, 9] = [0.000,0.000,0.000,0.000] ; 1997 CORY[0:3,12, 9] = [0.000,0.000,0.000,0.000] ; 1998 CORY[0:3,13, 9] = [0.000,0.000,0.000,0.000] ; 1999 CORY[0:3,14, 9] = [1.000,1.000,1.000,1.000] ; 2000 ;CORY[0:3,15, 9] = [1.000,0.961,1.000,0.036] ; 2001 ;CORY[0:3,16, 9] = [1.000,0.879,1.000,0.036] ; 2002 ;CORY[0:3,17, 9] = [0.968,0.671,1.000,0.010] ; 2003 ;CORY[0:3,18, 9] = [0.950,0.576,1.000,0.030] ; 2004 ;CORY[0:3,19, 9] = [0.855,0.381,1.000,0.067] ; 2005 ;CORY[0:3,20, 9] = [0.772,0.245,1.000,0.100] ; 2006 CORY[0:3,15, 9] = [1.000,0.961,0.983,1.000] ; 2001 CORY[0:3,16, 9] = [1.000,0.879,0.911,0.271] ; 2002 CORY[0:3,17, 9] = [0.968,0.671,0.907,0.831] ; 2003 CORY[0:3,18, 9] = [0.950,0.576,0.955,1.843] ; 2004 CORY[0:3,19, 9] = [0.852,0.378,0.841,2.564] ; 2005 CORY[0:3,20, 9] = [0.772,0.245,0.767,1.000] ; 2006 CORY[0:3,21, 9] = [0.563,0.097,0.729,1.000] ; 2007 CORY[0:3,22, 9] = [0.315,0.021,0.672,1.000] ; 2008 CORY[0:3, 0,10] = [0.000,0.000,0.000,0.000] ; F16 CORY[0:3, 1,10] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2,10] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3,10] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4,10] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5,10] = [0.000,0.000,0.000,0.000] ; 1991 CORY[0:3, 6,10] = [0.000,0.000,0.000,0.000] ; 1992 CORY[0:3, 7,10] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8,10] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9,10] = [0.000,0.000,0.000,0.000] ; 1995 CORY[0:3,10,10] = [0.000,0.000,0.000,0.000] ; 1996 CORY[0:3,11,10] = [0.000,0.000,0.000,0.000] ; 1997 CORY[0:3,12,10] = [0.000,0.000,0.000,0.000] ; 1998 CORY[0:3,13,10] = [0.000,0.000,0.000,0.000] ; 1999 CORY[0:3,14,10] = [0.000,0.000,0.000,0.000] ; 2000 CORY[0:3,15,10] = [0.000,0.000,0.000,0.000] ; 2001 CORY[0:3,16,10] = [0.000,0.000,0.000,0.000] ; 2002 CORY[0:3,17,10] = [1.000,1.000,1.000,1.000] ; 2003 CORY[0:3,18,10] = [1.000,1.000,1.000,1.000] ; 2004 ;CORY[0:3,19,10] = [0.842,0.841,0.832,0.832] ; 2005 ;CORY[0:3,20,10] = [0.688,0.685,0.628,0.628] ; 2006 CORY[0:3,19,10] = [0.593,0.593,0.593,0.593] ; 2005 CORY[0:3,20,10] = [0.567,0.567,0.567,0.567] ; 2006 CORY[0:3,21,10] = [0.463,0.463,0.463,0.463] ; 2007 CORY[0:3,22,10] = [0.460,0.460,0.460,0.460] ; 2008 CORY[0:3, 0,11] = [0.000,0.000,0.000,0.000] ; F17 CORY[0:3, 1,11] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2,11] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3,11] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4,11] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5,11] = [0.000,0.000,0.000,0.000] ; 1991 CORY[0:3, 6,11] = [0.000,0.000,0.000,0.000] ; 1992 CORY[0:3, 7,11] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8,11] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9,11] = [0.000,0.000,0.000,0.000] ; 1995 CORY[0:3,10,11] = [0.000,0.000,0.000,0.000] ; 1996 CORY[0:3,11,11] = [0.000,0.000,0.000,0.000] ; 1997 CORY[0:3,12,11] = [0.000,0.000,0.000,0.000] ; 1998 CORY[0:3,13,11] = [0.000,0.000,0.000,0.000] ; 1999 CORY[0:3,14,11] = [0.000,0.000,0.000,0.000] ; 2000 CORY[0:3,15,11] = [0.000,0.000,0.000,0.000] ; 2001 CORY[0:3,16,11] = [0.000,0.000,0.000,0.000] ; 2002 CORY[0:3,17,11] = [0.000,0.000,0.000,0.000] ; 2003 CORY[0:3,18,11] = [0.000,0.000,0.000,0.000] ; 2004 CORY[0:3,19,11] = [0.000,0.000,0.000,0.000] ; 2005 CORY[0:3,20,11] = [1.000,1.000,1.000,1.000] ; 2006 CORY[0:3,21,11] = [1.000,1.000,1.000,1.000] ; 2007 CORY[0:3,22,11] = [0.996,0.996,0.963,0.963] ; 2008 CORY[0:3, 0,12] = [0.000,0.000,0.000,0.000] ; F18 CORY[0:3, 1,12] = [0.000,0.000,0.000,0.000] ; 1987 CORY[0:3, 2,12] = [0.000,0.000,0.000,0.000] ; 1988 CORY[0:3, 3,12] = [0.000,0.000,0.000,0.000] ; 1989 CORY[0:3, 4,12] = [0.000,0.000,0.000,0.000] ; 1990 CORY[0:3, 5,12] = [0.000,0.000,0.000,0.000] ; 1991 CORY[0:3, 6,12] = [0.000,0.000,0.000,0.000] ; 1992 CORY[0:3, 7,12] = [0.000,0.000,0.000,0.000] ; 1993 CORY[0:3, 8,12] = [0.000,0.000,0.000,0.000] ; 1994 CORY[0:3, 9,12] = [0.000,0.000,0.000,0.000] ; 1995 CORY[0:3,10,12] = [0.000,0.000,0.000,0.000] ; 1996 CORY[0:3,11,12] = [0.000,0.000,0.000,0.000] ; 1997 CORY[0:3,12,12] = [0.000,0.000,0.000,0.000] ; 1998 CORY[0:3,13,12] = [0.000,0.000,0.000,0.000] ; 1999 CORY[0:3,14,12] = [0.000,0.000,0.000,0.000] ; 2000 CORY[0:3,15,12] = [0.000,0.000,0.000,0.000] ; 2001 CORY[0:3,16,12] = [0.000,0.000,0.000,0.000] ; 2002 CORY[0:3,17,12] = [0.000,0.000,0.000,0.000] ; 2003 CORY[0:3,18,12] = [0.000,0.000,0.000,0.000] ; 2004 CORY[0:3,19,12] = [0.000,0.000,0.000,0.000] ; 2005 CORY[0:3,20,12] = [0.000,0.000,0.000,0.000] ; 2006 CORY[0:3,21,12] = [0.000,0.000,0.000,0.000] ; 2007 CORY[0:3,22,12] = [1.000,1.000,1.000,1.000] ; 2008 ; ; GENERATE 4 DIGIT YEAR IF CALLED WITH 2 DIGIT AND DETERMINE YEAR ; INDEX INTO CORY MATRIX IYR=fix(NYR) IF IYR LT 200 then IYR=IYR+1900 IF IYR LT 1950 then IYR=IYR+100 IYX = IYR - 1986 IF IYX GT 20 THEN IYX = 20 ; Using last available correction IF IYX LT 0 THEN BEGIN print,' GETJGF CALLED WITH NYR = ',NYR,': NO CORRECTION AVAILABLE' IYX=0 ENDIF ; ; IDENTIFY SATELLITE ISAT=-1 IF IDSAT GT 100 THEN BEGIN ; This subroutine was called with 4-Digit Designator FOR ii=0,11 do begin IF IDSAT EQ IDS[II] then ISAT=II+6 ENDFOR ENDIF ELSE ISAT=IDSAT IF ( ISAT LT 6 OR ISAT GT 18 ) THEN BEGIN print,' GETJGF CALLED WITH IDSAT = ',IDSAT,': NOT RECOGNIZED' STOP ENDIF ; ; TRANSFER CHANNEL ENERGIES AND UNCORRECTED GEOMETRIC FACTORS TO RETURN VECTORS JJ=2 FOR II=0,22 DO BEGIN IF II LE 21 THEN BEGIN EE[II]=EGY[II] EI[II]=EGY[II] ENDIF IF II LE 2 THEN BEGIN GME[II]=1. GMI[II]=1. ENDIF ELSE BEGIN IF II NE 13 THEN BEGIN JJ=JJ+1 GME[JJ]=GG[II-3,0,ISAT-6] GMI[JJ]=GG[II-3,1,ISAT-6] ENDIF ENDELSE ENDFOR G11=GG[10,0,ISAT-6] FOR II=0,1 DO BEGIN FOR JJ=0,5 DO BEGIN ZF[JJ,II] = 1.0 IF ISAT GT 15 THEN ZF(JJ,II)=ZFAC[JJ,II,ISAT-16] ENDFOR ENDFOR ; ; APPLY CORY CORRECTION IF APPLICABLE FEH=1.0 FEL=1.0 FIH=1.0 FIL=1.0 CF[0]=1.0 CF[1]=1.0 CF[2]=1.0 CF[3]=1.0 IF (FIX(IFLG) MOD 2) NE 0 THEN BEGIN FEH= CORY(0,IYX,ISAT-6) < 1.0 ;The "less than" sign (<) is the IDL minimum operator. FEL= CORY(1,IYX,ISAT-6) < 1.0 FIH= CORY(2,IYX,ISAT-6) < 1.0 FIL= CORY(3,IYX,ISAT-6) < 1.0 CF[0]=FEH CF[1]=FEL CF[2]=FIH CF[3]=FIL IFG=0 IF FEH LT 1. AND FEH GT .001 THEN BEGIN FOR II=3,12 DO BEGIN IF( ISAT GE 16) THEN BEGIN FF=1.-(1.-FEH)*(II-4)/18.0 GME[II]=GME[II]*(1.-(1.-FEH)*(II-4)/18.0) GE1=GME[3] ENDIF ELSE BEGIN GME(II)=GME(II)*(.75+FEH/4.0) ENDELSE ENDFOR IFG=1 ENDIF IF FEL LT 1. AND FEL GT .001 THEN BEGIN FOR II=13,21 DO BEGIN IF( ISAT GE 16) THEN BEGIN GME(II)=GME(II)*(1.-(1.-FEL)*(II-4)/18.0) ENDIF ELSE BEGIN GME(II)=GME(II)*FEL ENDELSE ENDFOR IFG=1 ENDIF IF FIH LT 1. AND FIH GT .001 THEN BEGIN FOR II=3,12 DO BEGIN IF( ISAT GE 16) THEN BEGIN GMI(II)=GMI(II)*(1.-(1.-FIH)*(II-4)/18.0) ENDIF ELSE BEGIN GMI(II)=GMI(II)*(.75+FIH/4.0) ENDELSE ENDFOR IFG=1 ENDIF IF FIL LT 1. AND FIL GT .001 THEN BEGIN FOR II=13,21 DO BEGIN IF( ISAT GE 16) THEN BEGIN GMI(II)=GMI(II)*(1.-(1.-FIL)*(II-4)/18.0) ENDIF ELSE BEGIN GMI(II)=GMI(II)*FIL ENDELSE ENDFOR IFG=1 ENDIF ENDIF if (IFLG AND '80'X) NE 0 then begin GGE1=GME[3] XE1=XEIFC[ISAT-6] CF[0]=CF[0]*XEIFC[ISAT-6] CF[1]=CF[1]*XEIFC[ISAT-6] CF[2]=CF[2]*XIIFC[ISAT-6] CF[3]=CF[3]*XIIFC[ISAT-6] for II=3,21 do begin GME[II]=GME[II]*XEIFC[ISAT-6] GMI[II]=GMI[II]*XIIFC[ISAT-6] endfor endif RETURN END ;------------------------------------------------------------------------------ ; PRO raw2ephem - convert ephemeris data from packed integers to unpacked ; values (integer or float as necessary) ;------------------------------------------------------------------------------ PRO raw2ephem,bdata,jday,hour,minute,isecond,iyear,glat,glon,alt,gglat110,gglon110,$ ; 10/08 Emery: add cglon110 ; cglat110,mlthr,mltmin,mltsec cglat110,cglon110,mlthr,mltmin,mltsec jday = bdata(0) hour = bdata(1) minute = bdata(2) isecond = bdata(3) ; always = zero iyear = 1950 + bdata(4) glat = bdata(5) if glat gt 1800 then glat = float(glat-4995)/10. else glat = float(glat-900)/10. glon = float(bdata(6))/10. alt = bdata(7) gglat110 = float(bdata(8)) if gglat110 gt 1800 then gglat110 = float(gglat110-4995)/10. else $ gglat110 = float(gglat110-900)/10. gglon110 = float(bdata(9))/10. cglat110 = bdata(10) if cglat110 gt 1800 then cglat110 = float(cglat110-4995)/10. else $ cglat110 = float(cglat110-900)/10. cglon110 = float(bdata(11))/10. mlthr = bdata(12) mltmin = bdata(13) mltsec = bdata(14) end ;------------------------------------------------------------------------------ ; convert2counts - convert one minute of raw data to counts ; ; Parameters: ; bdata - array of raw data ; raw2count - array of values for converting tm data to counts ; counts - array of raw data converted to counts ; isatnumber - DMSP flight number, e.g. for F15, isatnumber=15 ; mode_flag - array of mode flags for 1 minute for SSJ5, 0=ModeA, 1=ModeB ; modeb_cycle - array of mode cycles for 1 min for SSJ5 in Mode B ; modeb_zone - value of Mode B zone for J5; set to zero for Mode A ;------------------------------------------------------------------------------ PRO convert2counts, bdata, raw2count, counts, isatnumber, mode_flag, modeb_cycle,$ modeb_zone ; Process 1 minute of data isectype = bdata[2595] ; Set to 1 for files processed after Sept 97 to indicate ; that times are stored as milliseconds instead of seconds for i=0,59 do begin ; Get the time tag for the second of data ; MilliSeconds of the Minute comes before 1 second of data isec = long(bdata[43*i+17]) if isectype ne 0 then begin if isec lt 0 then isec = 32768 + (isec and 32767) fsec = float(isec)/1000. isec = fix(fsec) endif ; Missing seconds of data are stored at end of 1 minute of data if i gt 0 and isec le 0 then goto,DONECVT ;Now convert Raw Data to Counts for j=0,9 do begin for k=0,3 do begin ; Done in groups of 4 because data stored as 4,3,2,1,8,7,6,5,etc. kk = bdata[43*i+18+4*j+3-k] if isec ge 0 and isec lt 60 then counts[4*j+k,isec] = raw2count[kk] ; electrons then ions endfor endfor ; For SSJ5, Decode Channel #E11 and #I11 if isatnumber gt 15 then begin ; Decode word E11 to find status ; MSB LSB ; 8 7 6 5 4 3 2 1 0 ; Z Z Z PA PA PA MD1 MD2 BD ; ; where ; Z's represents a 3 bit value of the Zone Address which is always zero in Mode A ; PA's represent a 3 bit value for the Parameter Address (whatever that means) ; MD1 is a 1 bit value for Mode 1 Select (whatever that means) ; MD2 is a 1 bit value for Mode 2 Select (whatever that means) ; BD is a 1 bit value for Boundary Detect s1 = bdata[43*i+18+9] ; Raw Data Word 11 = E11 za = fix(s1/64) ; Zone Address modeb_zone = za pa = fix((s1 - za*64)/8) ; Parameter Address md1 = fix((s1 - za*64 -pa*8)/4) ; Mode 1 Select md2 = fix((s1 - za*64 - pa*8 -md1*4)/2) ; Mode 2 Select bd = fix(s1 mod 2) ; Boundary Detect ; check Channel #E11 for whether instrument is in Mode A or B if za eq 0 then mode_flag[i]=0 else begin mode_flag[i] = 1 ; Mode B modeb_cycle[i] = za ; Sector counts are taken from endelse ; Decode word I11 to find status s2 = bdata[43*i+18+29] ; Raw Data Word 31 = I11 id = fix(s2/256) if pa eq 0 then begin wd = fix(s2/128)-256*id tp = fix(s2/64)-128*wd-256*id ebc = fix(s2/8)-64*tp-128*wd-256*id ibc = fix(s2 mod 8) endif else begin pv = fix(s2 mod 256) endelse endif endfor DONECVT: end ;------------------------------------------------------------------------------ ; PRO convert2flux - convert counts to differential number fluxes and store in ; array used for plotting ;------------------------------------------------------------------------------ ; PRO convert2flux, counts, num_min_saved_for_plot, $ e_fluxes, i_fluxes, e_geomf, i_geomf, isatnumber ;------------------------------------------------------------------------------ ; Parameters ; counts - one minute of counts from J4/J5 ; num_min_saved_for_plot - number of minutes of data that have been already ; saved for plotting ; e_fluxes - differential number fluxes of electron for plot interval ; i_fluxes - differential number fluxes of ions for plot interval ; e_geomf - geometric factors for 20 electron channels (30 KeV to 30 ev) ; i_geomf - geometric factors for 20 ion channels (30 KeV to 30 ev) ;------------------------------------------------------------------------------ nfluxes = fltarr(40,60) ; Convert from counts to differential fluxes for i=0,59 do begin if isatnumber le 15 then begin ; SSJ4 for j=0,19 do begin nfluxes[j,i] = counts[j,i]/(98.*e_geomf[j]) nfluxes[j+20,i] = counts[j+20,i]/(98.*i_geomf[j]) endfor endif else begin ; SSJ5 Treat Mode B as if it were mode A. If the flux ; is isotropic, a single geometric factor will give a small variation in ; flux over six second because the Mode B geometric factors vary slightly ; from sector to sector. for j=0,8 do begin nfluxes[j,i] = counts[j,i]/(50.*e_geomf[j]) nfluxes[j+11,i] = counts[j+11,i]/ (50.*e_geomf[j+11]) nfluxes[j+20,i] = counts[j+20,i]/(50.*i_geomf[j]) nfluxes[j+31,i] = counts[j+31,i]/ (50.*i_geomf[j+11]) endfor ; For SSJ5 only Channel 10 (and 30) are meaningful. Channels 11 (and 31) are ; status words. nfluxes[9,i] = counts[9,i]/(50.*e_geomf[9]) nfluxes[29,i] = counts[29,i]/(50.*i_geomf[9]) endelse endfor ; store the fluxes for i=0,59 do begin for j=0,19 do begin e_fluxes(60*num_min_saved_for_plot+i,j) = nfluxes(j,i) i_fluxes(60*num_min_saved_for_plot+i,j) = nfluxes(20+j,i) endfor endfor end PRO jbackgroundremoval,jcntin,jcntout,isat,abkg,itype ;----------------------------------------------------------------------------- ; pro jbackgroundremoval - procedure for identifying and removing background ; count rates from J4/J5 data. Since there is no back- ; ground channel, assume data is due to background when ; counts in most channels are approximately equal ; ; Parameters: ; jcntin - 20 channels of counts without background removal ; fltarr(20) ; 20 channels (elec or ions) ; jcntout - 20 channels of counts with background removal ; fltarr(20) ; 20 channels (elec or ions) ; isat - DMSP spacecraft number ; abkg - background count rate from period second of data ; fltarr(2) ; 2 species (elec, ions) ; itype - index for type of particle, 0=electrons, 1=ions ; ; Comments: ; Created January 2004 by Ernie Holman, AFRL/VSBX and Boston College. ; Looks for a nearly constant count rate in the higher energy channels of ; the electron or ion spectrum. If such a nearly constant count rate is ; found, then that rate plus one sigma of the rate is subtracted from all ; counts in the given spectrum. ;----------------------------------------------------------------------------- SX=0. SXX=0. NS=0 ; NJ = NUMBER OF CHANNELS USED FOR COMPUTING BACKGROUND ; FOR J4, USE 30KEV TO 1 KEV (CHANNELS 1 TO 10) ; FOR J5, USE 30KEV TO 100 EV (CHANNELS 1 TO 17) NJ=10 IF ISAT GT 15 THEN NJ=17 ; COMPUTE MEAN AND MEAN SQ OF COUNTS IN NJ CHANNELS for jj=0,nj-1 do begin if JJ NE 10 THEN begin SX = SX +jcntin[JJ] SXX = SXX+jcntin[JJ]*jcntin[JJ] NS=NS+1 ENDIF endfor IF NS EQ 0 then NS=1 AVG=SX/NS ; COMPUTE STD DEVIATION OF COUNTS FROM AVG SIGSQ = SXX/NS - (AVG^2) SIG = 0. SIGP = 0. IF SIGSQ GT 0 then SIG = SQRT(SIGSQ) ; COMPUTE PERCENTAGE DEVIATION IF AVG GT 1.E-6 then SIGP=SIG/AVG ; COMPUTE 1/AVG SIGQ=1. IF AVG GT 1.E-6 THEN SIGQ=1./SQRT(AVG) ; IF PERCENTAGE DEV IS LESS THAN 150 PERCENT OF POISSON ; STATISTICS VARIATION, THE COMPUTE BACKGROUND LEVEL AS ; AVG COUNT PLUS STD DEV, ; ELSE, USE VALUE FROM LAST ENTRY INTO THIS SUBROUTINE ; (ASSUMES PERSISTENT OF VALUES FROM ONE CALL TO THE ; NEXT) CUTS = 1.5*SIGQ IF SIGP LT CUTS THEN abkg[itype]=AVG+SIG ; DETERMINE HOW MANY VALUES WILL BE NEGATIVE IF USING ; INITIAL BACKGROUND VALUE NNEG=0 SNEG=0. ANEG=0. FOR JJ=0,19 DO BEGIN IF JJ NE 10 THEN BEGIN CCTMP = JCNTIN[JJ]-abkg[itype] IF CCTMP LT 0. THEN BEGIN NNEG=NNEG+1 SNEG=SNEG+CCTMP ENDIF ENDIF ENDFOR ;COMPUTER AVG OF NEG VALUES USING INITIAL BACKGROUND IF NNEG NE 0 THEN ANEG=SNEG/NNEG ;IF USING BACKGROUND FROM LAST SECOND INSTEAD OF ONE ;COMPUTED FOR THIS SECOND, DECREASE BACKGROUND BY HALF ;ON AVG NEG VALUE WHEN USING UN-MODIFIED BACKGROUND IF SIGP GE CUTS THEN abkg[itype]=abkg[itype]+ANEG/2. ; REMOVE BACKGROUND. NO COUNT CAN BE LESS THAN ZERO FOR JJ=0,19 DO BEGIN IF ISAT LE 15 OR JJ NE 10 THEN BEGIN jcntout[JJ]=jcntin[JJ] - abkg[itype] if jcntout[JJ] lt 0. then jcntout[JJ]=0. ENDIF ENDFOR ; CHECK FOR SINGLE CHANNEL WITH COUNTS SURROUNDED BY ; CHANNELS WITH ZERO COUNTS. IF FOUND, SET CHANNEL TO ; ZERO. THIS GETS RID OF RANDOM BAD DATA POINTS. for JJ=1,18 do begin IF JJ NE 10 THEN begin JM=JJ-1 JP=JJ+1 IF JJ EQ 9 then JP=11 IF JJ EQ 11 then JM=9 IF jcntout[JM] LT 0.1 AND jcntout[JP] LT 0.1 then jcntout[JJ] = 0. ENDIF endfor IF jcntout[18] LT 0.1 then jcntout[19]=0. IF jcntout[1] LT 0.1 then jcntout[0]=0. end ; End of procedures related to particle.pro ;------------------------------------------------------------------------------ ; PRO open_output_device - open the device upon which the graphics will be drawn ;------------------------------------------------------------------------------ PRO open_output_device, outtype, psfile,datafile_or_keyboard_input,$ ps_filename_problem print,"Enter destination of plot." STARTSET: if datafile_or_keyboard_input eq 0 then begin print,' Enter device no. according to this scheme: ' print,' 1) Postscript file (unix)' print,' 2) Hewlett-Packard Printer file (unix)' print,' 3) X terminal (unix)' print,' 4) Z buffer --> *.png file (unix)' print,' -1) Postscript file (win/vms)' print,' -2) Hewlett-Packard Printer file (win/vms)' print,' -3) MS Windows plot (win/vms)' print,' -4) Z buffer --> *.png file (win/vms)' read, outtype endif else readf,1,outtype print,outtype !P.FONT =-1 !P.THICK = 1.0 if outtype eq -3 then begin set_plot,'WIN' device,decomposed=0 ; Use 8-bit graphics on a 16-bit or 24-bit graphics system goto, SETDONE endif if outtype eq 3 then begin set_plot, 'X' ; X-windows plot device,pseudo_color=8 goto, SETDONE endif if abs(outtype) eq 1 then begin if datafile_or_keyboard_input eq 0 then begin print,'Indicate the type of post-script file desired by entering:' print,' 0 for standard post-script file,' print,' 1 for encapsulated post-script file' read,ps_encapsulated_option endif else begin readf,1,ps_encapsulated_option endelse if ps_encapsulated_option lt 0 or ps_encapsulated_option gt 1 then $ ps_encapsulated_option=0 psfile = " " GET_PS_FILE_NAME: if datafile_or_keyboard_input eq 0 then begin if ps_encapsulated_option eq 0 then $ print,'Enter filename for the post-script output file. (include .ps extension)' else $ print,'Enter filename for the encapsulated post-scipt output file. (include .eps extension)' read,format='(a)',psfile endif else readf,1,format='(a)',psfile print,psfile on_ioerror,PSFILE_NOT_FOUND get_lun,u_ps openr,u_ps,psfile ; Can only get to the next line if psfile already exists close,u_ps print,'*** Post Script File Name Is Already Being Used ***' if datafile_or_keyboard_input ne 0 then begin print,'IDL run is aborting to avoid over-writing the file.' free_lun,u_ps OUTTYPE = -99 goto,SETDONE endif print,'Enter an integer to indicate what to do next' print,' 1 = Over-write file' print,' 2 = Get new file name' print,' 3 = Abort IDL run' read,ps_filename_problem if ps_filename_problem eq 2 then begin free_lun,u_ps GOTO,GET_PS_FILE_NAME endif if ps_filename_problem eq 3 then begin free_lun,u_ps outtype = -99 goto,SETDONE endif PSFILE_NOT_FOUND: free_lun,u_ps ; A satisfactory postscript filename has been found. goto, SETDONE endif if abs(outtype) eq 2 then begin print,"Enter the path/name of the output file" fn = " " read,format='(a)',fn print,fn set_plot, 'PCL'; HP PCL file !P.THICK = 3.0 print,"PCL" device,/landscape,filename=fn goto, SETDONE endif if abs(outtype) eq 4 then begin print,"Z-buffer" set_plot,'Z' ; Z-buffer device,set_resolution=[800,600] goto, SETDONE endif ; Error - unknown file format outtype = -99 SETDONE: return end ;----------------------------------------------------------------------------- ; ; procedure get_file_names ; ; Reads in particle data file name, driftmeter file name and ; the postscript file name. ; ;----------------------------------------------------------------------------- PRO get_file_names, partfile$, removebackground, driffile$, ssmffile$, smfile$,$ psfile, datafile_or_keyboard_input, output_type, dm_or_ssm_option, $ inputerror, ps_encapsulated_option inputerror = 0 datafile_or_keyboard_input = 0 print, 'Input information from keyboard or data file (k or d) ?' dum$ = '' read, dum$ if strmid(dum$,0,1) eq 'd' or strmid(dum$,0,1) eq 'D' then $ datafile_or_keyboard_input=1 $ else datafile_or_keyboard_input=0 if datafile_or_keyboard_input eq 1 then begin print, 'Enter filename containing input stream parameters :' fileinput = '' on_ioerror,BADINPUTFILE read, fileinput openr,1,fileinput on_ioerror,NULL goto,STARTREADINGINPUT BADINPUTFILE: on_ioerror,NULL print,'Filename for inputs not found. End of Job.' inputerror = -1 goto,GETFILEDONE endif ; ; Get the index number for the Output Device (variable name is 'output_type') STARTREADINGINPUT: output_type=-99 psfile = '' ps_filename_problem = 0 open_output_device, output_type,psfile,datafile_or_keyboard_input,ps_filename_problem if output_type le -99 then goto,GETFILEDONE ; Start getting the name of the data files partfile$ = '' driffile$ = 'null' ssmffile$ = 'null' smfile$ ='null' ; print, 'Enter path/Filename for Particle data (null for no data):' if datafile_or_keyboard_input eq 0 then read, partfile$ else readf,1,partfile$ print, 'The name of the file selected is : ',partfile$ on_ioerror,WRONGJ4FILE GET_LUN,UJ4 OPENJ4FILE: OPENR,UJ4,partfile$ close,uj4 goto,J4FILEOK WRONGJ4FILE: PRINT,"Could not open J4/J5 data file: ",partfile$ if datafile_or_keyboard_input eq 1 then begin free_lun,uj4 inputerror = -1 print,'IDL job aborting' close,/all goto,GETFILEDONE endif PRINT,'Re-Enter path/Filename for Particle data (null for no data):' read,partfile$ if partfile$ ne "null" and partfile$ ne "NULL" then goto,OPENJ4FILE J4FILEOK: on_ioerror,NULL free_lun,uj4 print,'Indicate whether the J4/J5 background should be removed.' print,' Enter 0 for NO, and any value NE 0 for YES' if datafile_or_keyboard_input eq 0 then read, removebackground else $ readf, 1, removebackground print,removebackground right$ = 'n' print, 'Indicate choice of driftmeter, Ni and/or magnetometer data.' print,' 1) Driftmeter data in lower box.' print,' 2) Magnetometer data in lower box.' print,' 3) Lower box empty.' print,' 4) Driftmeter data in lower box. SSM Delta-B & SM Ni in upper boxes' print,' 5) Driftmeter data in lower box. SSM Delta-B in larger upper box' if datafile_or_keyboard_input eq 0 then read, dm_or_ssm_option else readf, 1, dm_or_ssm_option if dm_or_ssm_option ne 1 and dm_or_ssm_option ne 2 and dm_or_ssm_option ne 4 $ and dm_or_ssm_option ne 5 then dm_or_ssm_option = 3 print, dm_or_ssm_option ; ; Get Driftmeter file name if dm_or_ssm_option eq 1 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin print,'Enter path/filename for Driftmeter data (null for no data):' if datafile_or_keyboard_input eq 0 then read,driffile$ else readf,1,driffile$ print,driffile$ print, 'The name of the IDM file selected is : ',driffile$ if driffile$ eq "null" or driffile$ eq "NULL" then goto,I2FILEOK on_ioerror,WRONGI2FILE GET_LUN,UI2 OPENI2FILE: OPENR,UI2,driffile$ close,UI2 goto,I2FILEOK WRONGI2FILE: PRINT,"Could not open DM file: ",driffile$ if datafile_or_keyboard_input eq 0 then begin PRINT,"Re-Enter Filename for Driftmeter data (null for no data):' read,driffile$ if driffile$ ne "null" and driffile$ ne "NULL" then goto,OPENI2FILE endif else begin inputerror = -1 print,'IDL job aborting' close,/all goto,GETFILEDONE endelse I2FILEOK: on_ioerror,NULL free_lun,UI2 endif ; ; Get SSM (magnetometer) file name if dm_or_ssm_option eq 2 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin print,'Enter path/filename for Magnetometer data (null for no data):' if datafile_or_keyboard_input eq 0 then read,ssmffile$ else readf,1 ,ssmffile$ print,'The name of the SSM file is: ',ssmffile$ if ssmffile$ eq "null" or ssmffile$ eq "NULL" then goto,MFILEOK on_ioerror,WRONGMFILE GET_LUN,UM OPENMFILE: OPENR,UM,ssmffile$ close,UM goto,MFILEOK WRONGMFILE: PRINT,"Could not open SSM file: ",ssmffile$ if datafile_or_keyboard_input eq 1 then begin inputerror = -1 print,'IDL job aborting' close,/all goto,GETFILEDONE endif PRINT,"Re-Enter Filename for Magnetometer data (null for no data):' read,ssmffile$ if ssmffile$ ne "null" and ssmffile$ ne "NULL" then goto,OPENMFILE MFILEOK: on_ioerror,NULL free_lun,UM endif ; ; Get SM (SSIES Total Ion Trap) file name if dm_or_ssm_option eq 4 then begin print,'Enter path/filename for SSIES SM data (null for no data):' if datafile_or_keyboard_input eq 0 then read,smfile$ else readf,1 ,smfile$ print,'The name of the SSIES/SM file is: ',smfile$ if smfile$ eq "null" or smfile$ eq "NULL" then goto,SMFILEOK on_ioerror,WRONGSMFILE GET_LUN,USM OPENSMFILE: OPENR,USM,smfile$ close,USM goto,SMFILEOK WRONGSMFILE: PRINT,"Could not open SM file: ",smfile$ if datafile_or_keyboard_input eq 1 then begin inputerror = -1 print,'IDL job aborting' close,/all goto,GETFILEDONE endif PRINT,"Re-Enter Filename for SM data (null for no data):' read,smfile$ if smfile$ ne "null" and smfile$ ne "NULL" then goto,OPENSMFILE SMFILEOK: on_ioerror,NULL free_lun,USM endif if dm_or_ssm_option eq 3 then begin ; Lower box is empty driffile$ = 'null' ssmffile$ = 'null' smfile$ = 'null' endif ; GETFILEDONE: RETURN END ;----------------------------------------------------------------------------- ; ; procedure get_time_int ; ; The procedure simply reads in all of the other inputs for the run. ; ; BEGTIME - array of integers for beginning time (see above) ; ENDTIME - time in seconds for the plot to end ; LENGTH - maximum time in minutes for any plot ; SATEL - satellite number ; COR - remove corotation = 0.0, not = 1.0 ; YCOL - color = 1.0, b/w = 0.0 ; YFIL - filter ave energy and flux = 0.0, no = 1.0 ; MAXVEL - maximum velocity to plot (set's range) ; MINLAT - minimum latitude to plot ; set to 0.0 if NODM = 1 ; NSQ - if 1, northern hemisphere ; - if -1.0, southern hemisphere ; NODM - if 0, plot dm in lower box ; - if 1, value of drffile = 'null' ; NOSSM - if 0, plot ssm in lower box ; - if 1, value of ssmfile = 'null' ; VOPT - if 0, do not plot the vertical driftmeter ; - if 1, yes, plot the vertical driftmeter ; datafile_or_keyboard_input - if 0, read input from keyboard ; - if 1, readf input from a prepared file ; portrait_or_land_option - Landscape (0) portrait (1) ; ; DISWITCH IF 0, DISPLAY THE SPECTOGRAM AS THE NUMBER FLUX ; IF 1, DISPLAY THE SPECTOGRAM AS THE ENERGY FLUX ; ; INTSW IF 0, DISPLAY THE TOP B/W PLOT AS INTEGRATED NUMBER FLUX ; IF 1, DISPLAY THE TOP B/W PLOT AS INTEGRATED ENERGY FLUX ; ; IMUL MULTIPLIER (factor of 1,2 or 3 for spacing the tick marks) ; ; More comments to help with understanding the parameters: ; begtime(0) = year ; begtime(1) = month ; begtime(2) = day of month ; begtime(3) = beginning hour ; begtime(4) = beginning minute ; begtime(5) = julian day ;----------------------------------------------------------------------------- PRO get_time_int, dayofmon, begtime, endtime, length, satel,cor, ycol, yfil, $ maxvel, max_nt, minlat, nsq, nodm, nossm, vopt,$ datafile_or_keyboard_input, portrait_or_land_option, diswitch, flipi,$ vcount, vlabc, vstr, vlabx, vlines, intsw, imul, output_type if abs(output_type) eq 1 then begin portrait_or_land_option = '' print, ' ' print, 'Would you like landscape (l) or portrait (p) ?) if datafile_or_keyboard_input eq 0 then read, portrait_or_land_option else readf,1,portrait_or_land_option print, portrait_or_land_option if (strmid(portrait_or_land_option,0,1) eq 'p') or (strmid(portrait_or_land_option,0,1) eq 'P') then portrait_or_land_option = 1 $ else portrait_or_land_option = 0 endif else portrait_or_land_option = 0 ans = '' print, ' ' print, ' Would you like to display the spectogram as' print, ' differential (E)nergy Flux (ans: E) or ' print, ' differential (N)umber flux (ans: N) ?' if datafile_or_keyboard_input eq 0 then read, ans else readf,1, ans print,ans diswitch = 0 if (strmid(ans,0,1) eq 'E') or (strmid(ans,0,1) eq 'e') then diswitch = 1 ans = '' print, ' ' print, ' Would you like to display the TOP b/w section as ' print, ' integrated Energy Flux (ans: E) or ' print, ' integrated Number Flux (ans: N) ?' if datafile_or_keyboard_input eq 0 then read, ans else readf,1, ans print,ans intsw = 0 if (strmid(ans,0,1) eq 'E') or (strmid(ans,0,1) eq 'e') then intsw = 1 dum = 0 print, 'Enter Year of data to be plotted :' if datafile_or_keyboard_input eq 0 then read, dum else readf,1, dum print, dum ; Get the year from the user if (dum lt 100) then begin if dum gt 86 then begtime[0] = 1900+dum else begtime[0] = 2000+dum endif else begtime[0] = dum print, 'Enter Month of Data to be plotted :' if datafile_or_keyboard_input eq 0 then read, dum else readf,1, dum print,dum begtime[1] = dum print, 'Enter Day of the month to be plotted :' if datafile_or_keyboard_input eq 0 then read, dum else readf,1, dum print,dum begtime[2] = dum print, 'Enter Time of day to begin plotting (enter HHMM, HHMM:00 assumed ) :' if datafile_or_keyboard_input eq 0 then read, dum else readf,1, dum print,dum begtime[3] = fix(dum/100) begtime[4] = dum-(begtime(3)*100) if fix(float(begtime[0])/4.0) ne fix(float(begtime[0]-1)/4.0) then $ dayofmon(1) = 29 begtime[5] = 0 for dum = 0, begtime[1]-2 do begtime[5] = begtime[5] + dayofmon(dum) begtime[5] = begtime[5] + begtime[2] dum = 0 print, 'Enter Time of day to END plotting (enter HHMM, HHMM:59 assumed ) :' if datafile_or_keyboard_input eq 0 then read, dum else readf,1, dum print,dum enddum = fix(dum/100) endtime = float(enddum)*3600.0 + float(dum-(enddum*100))*60.0 print, ' ' print, 'Specify the length of time for each page (max 30 min):' if datafile_or_keyboard_input eq 0 then read, length else readf,1, length print,length if length gt 30. then begin print,'Selected length exceeds program limit of 30 minutes.' print,'Length is being re-set to 30 minutes/plot' length = 30. endif ; Comment: The SunOS 4.1 VERSION always forces a tic mark on the right ; side of the plot. That means there cannot be plot lengths of 20 or ; less that cannot be divided by 2 or plot lengths of 21-30 that cannot ; be divided by 3. So this code is code from procedure tick_stuff that ; has been modified and moved here to make the minutes be displayed ; correctly. if length le 11.0 then begin nlen = length imul = 1 endif else begin if length le 20.0 then begin nticks = fix(length/2.0 + 0.6) imul = 2 nlen = imul * nticks endif else begin nticks = fix(length/3.0 + 0.7) imul = 3 nlen = imul * nticks endelse endelse length = float(nlen) ; ; if nodm eq 0 then dum = -1.0 else dum = 0.0 dum = -1.0 while (dum gt 80.0) or (dum lt 0.0) do begin print, ' ' print, 'Enter minimum latitude to plot (80.0 Max) :' if datafile_or_keyboard_input eq 0 then read, dum else readf, 1, dum print, dum endwhile minlat = dum nsq = 0.0 dum = 0 while (dum gt 3) or (dum lt 1) do begin if datafile_or_keyboard_input eq 0 then begin print, ' ' print, '1. Northern Hemisphere only' print, '2. Southern Hemisphere only' print, '3. Both Hemispheres' print, 'Enter choice:' read, dum endif else readf,1, dum endwhile ; The value of dum is used to set nsq if dum eq 1 then nsq = 1.0 if dum eq 2 then nsq = -1.0 cor = 1.0 vopt = 0 ; This switch, nodm, gets set in the M A I N program according to whether ; driffile is null or a valid filename nodm = 0, nossm = 1 ; ; Another switch, nossm, also gets set in the MAIN program according to ; whether ssmffile is null or a valid file name, nossm = 0, nodm = 1 if nodm eq 0.0 then begin que$='' if datafile_or_keyboard_input eq 0 then begin print, 'Corotation to be removed from the Horizontal Velocity (y/n) ?' read, que$ endif else readf,1, que$ if (strmid(que$,0,1) eq 'y') or (strmid(que$,0,1) eq 'Y') then cor = 0.0 que$='' if datafile_or_keyboard_input eq 0 then begin print, 'Would you like the vertical driftmeter data plotted (y/n)?' read, que$ endif else readf,1, que$ if (strmid(que$,0,1) eq 'y') or (strmid(que$,0,1) eq 'Y') then vopt = 1 if datafile_or_keyboard_input eq 0 then begin print, ' ' print, 'Enter maximum Velocity for Drift Meter (3000.0 Max) :' read, maxvel endif else readf,1, maxvel endif ; if nossm eq 0 then begin if datafile_or_keyboard_input eq 0 then begin print, 'Enter maximum nana-tesla magnetometer measurement:' read, max_nt endif else readf,1,max_nt endif if datafile_or_keyboard_input eq 0 then begin print, 'Enter Satellite number (6-17) :' read, dum endif else readf,1, dum if dum lt 6 then dum=6 if dum gt 17 then dum = 17 satel = fix(dum) ycol = 1.0 que$='' if datafile_or_keyboard_input eq 0 then begin print, 'Would you like the plots to be in color (y/n) ?' read, que$ endif else readf,1, que$ if (strmid(que$,0,1) eq 'n') or (strmid(que$,0,1) eq 'N') then ycol = 0.0 yfil = 0 que$='' if datafile_or_keyboard_input eq 0 then begin print, 'Would you like the Energy Flux and Average Energy ' print, 'plots to be filtered (y/n) ?' read, que$ endif else readf,1, que$ if (strmid(que$,0,1) eq 'n') or (strmid(que$,0,1) eq 'N') then yfil = 1 flipi = 0 que$='' if datafile_or_keyboard_input eq 0 then begin print, 'Would you like to have the ion spectrogram upside down ? ' read, que$ endif else readf,1, que$ if (strmid(que$,0,1) eq 'y') or $ (strmid(que$,0,1) eq 'Y') then flipi = 1 lab = 0 vcount = 0 ; Number of vertical lines vlabc = 0 ; Number of labels vlines = fltarr(1) ; X value for vertical lines vstr = strarr(1) ; String for each label. Array size expands as user enters labels. vlabx = fltarr(1) ; X value for each label. que$='' if datafile_or_keyboard_input eq 0 then begin print, 'Would you like to have lines and/or labels on your plot? (y/n) ' read, que$ endif else readf,1, que$ if (strmid(que$,0,1) eq 'y') or $ (strmid(que$,0,1) eq 'Y') then lab = 1 if lab eq 1 then begin ; User asked for Lines and/or Labels. Get details. done = 0 while done eq 0 do begin que$='' if datafile_or_keyboard_input eq 0 then begin print, ' ' print, 'Enter time in HHMMSS to put vertical line and/or label,' print, 'or press enter to exit the label/line input :' read, que$ endif else readf,1, que$ if (strlen(que$) eq 0) or $ (min(byte(strmid(que$,0,2))) lt 48) or $ (max(byte(strmid(que$,0,2))) gt 57) then begin print, 'Exiting line/label section' done = 1 endif else begin hour = float(strmid(que$,0,2)) if strlen(strmid(que$,2,2)) gt 0 then $ min = float(strmid(que$,2,2)) $ else min = 0.0 if strlen(strmid(que$,4,2)) gt 0 then $ sec = float(strmid(que$,4,2)) $ else sec = 0.0 x = hour*3600.0+min*60.0+sec if (x lt 0.0) or (x gt 86400.0) then begin print, 'This is an invalid time.' print, 'Simply press < Return > to exit.' endif else begin que$='' if datafile_or_keyboard_input eq 0 then begin print, ' ' print, 'Would you like a Veritical Line at this time ?' read, que$ endif else readf,1, que$ if (strmid(que$,0,1) eq 'y') or $ (strmid(que$,0,1) eq 'Y') then begin if vcount eq 0 then vlines(0)=x else vlines=[vlines,x] vcount = vcount + 1 endif que$='' if datafile_or_keyboard_input eq 0 then begin print, ' ' print, 'Enter String for this time (return for no string) :' read, que$ endif else readf,1, que$ if (strlen(que$) gt 0) then begin if vlabc eq 0 then begin vstr(0)=que$ vlabx(0)=x endif else begin vstr=[vstr,que$] vlabx=[vlabx,x] endelse vlabc = vlabc + 1 endif endelse endelse endwhile endif if datafile_or_keyboard_input eq 1 then close,1 RETURN END ;----------------------------------------------------------------------------- ; ; procedure tick_stuff ; ; This procedure takes the amount of time to be plotted (LENGTH) and ; determines how many tick marks to have (NTICKS), where they should ; go (TICKVAL), the labels (TICKLAB) and the amount of time, in minutes, ; between each tick mark (IMUL - needed when writing out the ephemeris ; data). ; ; Each tick mark must land on an exact minute, since the ephemeris ; data is on the minute. ;----------------------------------------------------------------------------- PRO tick_stuff, begtime, length, nticks, tickval, ticklab, bt, et, imul ; Figure out the amount of ticks needed and the length of time between ; Each tick mark: ; if length le 11.0 then begin nticks = fix(length) dt = 60.0 imul = 1 endif else begin if length le 20.0 then begin nticks = fix(length/2.0 ) dt = 120.0 imul = 2 endif else begin nticks = fix(length/3.0 ) dt = 180.0 imul = 3 endelse endelse tickval = fltarr(nticks+1) ticklab = strarr(nticks+1) ; Set the first time to the beginning of the plot: utsec = float(begtime(3))*3600.0 + float(begtime(4))*60.0 for i = 0, nticks do begin ; Make the time into a string and call it TICKLAB ticklab(i) = strcompress(string(utsec),/remove_all) j = 0 ; We know that the time is in exact seconds, so cut it off at the ; decimal point: while strmid(ticklab(i),j,1) ne '.' do j=j+1 ticklab(i) = strmid(ticklab(i),0,j) ; Set the position of the tick mark: tickval(i) = float(i) * dt ; Increment the clock: utsec = utsec + dt endfor ; Set first time and last time on the plot: ; This is where bt and et are established and bt is always 0.0. bt = 0.0 et = length*60.0 RETURN END ;----------------------------------------------------------------------------- ; procedure set_color_palette ; ; ; creates a 200 point color table if color is selected ; ;----------------------------------------------------------------------------- PRO set_color_palette, output_type, ycol red = fltarr(256) green = fltarr(256) blue = fltarr(256) ; if ycol gt 0 then begin ; color table for spectrogram for i=1,39 do begin red[i] = 120. + (i*90.)/39. green[i] = 0.0 blue[i] = 140.+ float(i * 100) / 39. endfor for i=40,79 do begin red[i] = 200. - float(i-40) * 200./39. green[i] = float(i - 40) * 225./39. blue[i] = 240. - float(i - 39) / 2. endfor for i=80,105 do begin ; red[i] = 0.0 ; green[i] = 200.0 ; blue[i] = 200.0 - float(i-79)*200.0/26.0 ; red[i] = 0.0 green[i] = 225.0 - float(i - 80) blue[i] = 220. - float(i - 79) * 210./26. endfor for i=106,132 do begin red[i] = float(i-105)*255.0/27.0 green[i] = 200.0 + float(i-105)*55.0/27.0 blue[i] = 0.0 endfor for i=133,159 do begin red[i] = 255.0 green[i] = 240.0 - float(i-133)*240.0/26.0 blue[i] = 0.0 endfor for i=160,200 do begin red[i] = 255.0 green[i] = float(i-160)*5.0 blue[i] = float(i-160)*5.0 endfor for i=201,255 do begin red[i] = i green[i] = i blue[i] = i endfor endif else begin ; greyscale for spectrogram for i=0,255 do begin red[i] = i green[i] = i blue[i] = i endfor endelse ; ; POSTSCRIPT mode output_type = +/- 1 ; Hewlett-Packard Printer output_type = +/- 2 ; if abs(output_type) le 2 then begin red[0] = 0. green[0] = 0.0 blue[0] = 0. red[1] = 0 green[1] = 0 blue[1] = 1 endif ; ; X- WINDOW MODE output_type = 3 ; MS Windows Mode output_type = -3 ; Z-Buffer Mode output_type = +4 or -4 if abs(output_type) ge 3 then begin ; white red[0] = 255. green[0]= 255. blue[0] = 255. ; black red[1] = 0. green[1]= 0. blue[1] = 0. endif ; tvlct, red,green,blue RETURN END ;----------------------------------------------------------------------------- ; ; procedure plot_electron ; ; simply sets up the ranges for the electron SSJ/4 data, then plots ; it. ; ;----------------------------------------------------------------------------- PRO plot_electron, ele2d, bt, et, nticks, tickval, ticklab, yticklab, change,$ plotparsw, dm_or_ssm_option, ycol ; ele2d - array of color indices to indicate flux levels. For particles, ; index goes from 0 to 200 (color) and 255 (white) for background. ; bt - beginning time ; et - end time ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait ; plotparsw - flag to indicate whether particle data was found for time ; interval (=1) or was not found (=0) ; dm_or_ssm_option - flag to indicate whether DM data (=1), SSM data (=2), ; neither (=3) or both (=5) or both plus SM Ni data (=4) will ; be plotted with the particle data ; YCOL - flag to indicate whether plot is in color or b&w: color=1, b/w=0 ; if dm_or_ssm_option eq 3 then begin ymax_plote = 0.715 ymin_plote = 0.420 endif else begin ymax_plote = 0.715 ymin_plote = 0.520 endelse ysize_ele = ymax_plote - ymin_plote - 0.001 plot, [bt,et], [31.0,30500.0], $ /ylog, xstyle = 1, $ xticks = nticks, $ xtickname = ticklab, $ xtickv = tickval, $ xticklen = -0.015, $ ystyle = 1, $ ytickn = yticklab, $ yticklen = -0.005, $ charsize = change, $ pos = [0.09,ymin_plote,0.88,ymax_plote], $ color=1, /nodata if ycol lt 1 then begin ; B&W - re-map 1 to 200 into 40 to 240 in steps of 5 ; This increases the contrast between levels. esize = size(ele2d) eisize = esize[1] ejsize = esize[2] for i=0,eisize-1 do begin for j=0,ejsize-1 do begin if ele2d[i,j] gt 0 and ele2d[i,j] lt 201 then begin k = ele2d[i,j] k = fix(float(k)/5) ele2d[i,j] = fix(40. + 5*k ) endif endfor endfor endif ; The electron spectogram bitmap is "stretched" to fit the area specified by ; Xsize and Ysize, only on Postscript and HP output devices; not on X-windows, ; or MS-Windows devices. if plotparsw ne 0 then tv, ele2d, 0.0905, ymin_plote+0.001, xsize=0.789,$ ysize=ysize_ele,/norm RETURN END ;----------------------------------------------------------------------------- ; ; procedure plot_ion ; ; simply sets up the ranges for the ion SSJ/4 data, then plots it. ; ;----------------------------------------------------------------------------- PRO plot_ion, ion2d, bt, et, nticks, tickval, ticklab, yticklab, $ change, flipi, plotparsw, dm_or_ssm_option,ycol ; ion2d - array of color indices to indicate flux levels. For particles, ; index goes from 0 to 200 (color) and 255 (white) for background. ; bt - beginning time ; et - end time ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait ; plotparsw - flag to indicate whether particle data was found for time ; interval (=1) or was not found (=0) ; dm_or_ssm_option - flag to indicate whether DM data (=1), SSM data (=2), ; neither (=3) or both (=5) or both plus the SM Ni data will ; be plotted with the particle data ; YCOL - flag to indicate whether plot is in color or b&w: color=1, b/w=0 ; ; ytickn = strarr(3) if flipi eq 0 then begin ybot = 31.0 ytop = 30500.0 ytickn(0:2) = yticklab(0:2) endif else begin ytop = 31.0 ybot = 30500.0 ytickn[0] = yticklab(2) ytickn[1] = yticklab(1) ytickn[2] = yticklab(0) endelse if dm_or_ssm_option eq 3 then begin ymax_ploti = 0.415 ymin_ploti = 0.120 endif else begin ymax_ploti = 0.515 ymin_ploti = 0.320 endelse ysize_ion = ymax_ploti - ymin_ploti - 0.001 plot, [bt,et], [ybot, ytop], $ /ylog, xstyle = 1, $ xticks = nticks, $ xtickname = ticklab, $ xtickv = tickval, $ xticklen = -0.015, $ ytickn = ytickn, $ yticklen = -0.005, $ yrange = [ybot, ytop], $ ystyle = 1, $ charsize = change, $ pos = [0.09,ymin_ploti,0.88,ymax_ploti], $ /noerase, color=1, $ /nodata if ycol lt 1 then begin ; B&W - re-map 1 to 200 into 40 to 240 in steps of 5 ; This increases the contrast between levels. isize = size(ion2d) iisize = isize[1] ijsize = isize[2] for i=0,iisize-1 do begin for j=0,ijsize-1 do begin if ion2d[i,j] gt 0 and ion2d[i,j] lt 201 then begin k = ion2d[i,j] k = fix(float(k)/5) ion2d[i,j] = fix(40. + 5*k ) endif endfor endfor endif if flipi eq 1 then begin xr = et-bt for i= 40, 90, 10 do begin ; create tick marks between 100 eV and 30 eV plots, [-0.005*xr/2.0,0.0],[i,i] plots, [et+0.005*xr/2.0,0.0],[i,i] endfor if plotparsw ne 0 then tv, ion2d, 0.0905,ymin_ploti+0.001,xsize=0.789,$ ysize=ysize_ion,/norm, /order endif else $ if plotparsw ne 0 then tv, ion2d, 0.0905,ymin_ploti+0.001,xsize=0.789,$ ysize=ysize_ion,/norm RETURN END ;----------------------------------------------------------------------------- ; ; procedure plot_dm ; ; sets up ranges for the driftmeter data, then plots it out. This is ; a little more difficult than the other data, since it has to be plotted ; out with the time, so any data gaps (since the time values are set ; to -1 in the code which read the data file) will cause problems. So time ; values are searched to find where there are negative and positive steps. ; Positive steps mean that there is really data, negative steps mean ; that the end of the data has arrived or there is a data gap, since ; time progresses only positively and -1 implies data missing. ; So positive steps are plotted and negative steps are ignored. ; ; After the data is plotted, the ephemeris data is displayed ; in the subroutine plot_extra. ; ;----------------------------------------------------------------------------- PRO plot_dm, dmdata, bt, et, nticks, tickval, ticklab, $ yb, ye, ytick, ytickval, yticklab, mlt, $ mlat, begtime, npts, cor, imul, nodm, ycol, $ change, vopt ; ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait tickfake=strarr(20) for i=0,19 do tickfake(i)=' ' plot, [bt,et], [yb,ye], $ xstyle = 1, $ xticks = nticks, $ xtickname = tickfake, $ xtickv = tickval, $ xticklen = -0.03, $ ystyle = 1, $ yticks = ytick, $ ytickname = yticklab, $ ytickv = ytickval, $ yticklen = -0.005, $ yminor = 5, $ charsize = change, $ pos = [0.09,0.12,0.88,0.315], $ /noerase, color=1, $ /nodata if ycol eq 1.0 then begin hcolor = 1 ; Plot in color vcolor = 107 endif else begin ; Plot in B&W hcolor = 1 vcolor = 1 endelse hem = where(mlat ne 0.0, count) if count ge 2 then begin x = (90.0-abs(mlat))*sin(mlt*!pi/12.0) hem = where(mlat(0:1) gt 0.0, count) if count gt 0 then x = -1.0*x if x(1)-x(0) gt 0.0 then begin xyouts, et-(et-bt)/100.0, yb*0.75, 'Sunward', alignment=1.0, $ charsize=change,color=1 xyouts, et-(et-bt)/100.0, -yb*0.75, 'Antisunward', alignment=1.0, $ charsize=change,color=1 endif else begin xyouts, et-(et-bt)/100.0, -yb*0.75, 'Sunward', alignment=1.0, $ charsize=change,color=1 xyouts, et-(et-bt)/100.0, yb*0.75, 'Antisunward', alignment=1.0, $ charsize=change,color=1 endelse endif if nodm eq 0 then begin ; IDM data does exit. Plot it utoff = begtime(3)*3600.0+begtime(4)*60.0 ; get timesteps: dmsize = size(dmdata) dmd2 = dmsize[2] dumdat = dmdata(0,1:dmd2-1) - dmdata(0,0:dmd2-2) ; where are they positve: nptsarr = where(dumdat gt 0.0,count) ; If all the data is good, decide where the end of the data is, ; grab it and plot it: last = nptsarr(n_elements(nptsarr)-1)+1 timdata = dmdata(0,0:last)-utoff verdata = dmdata(1,0:last) hordata = dmdata(2,0:last) timdum = timdata-float(fix(timdata)) srch = where(timdum eq 0.0, scount) if vopt eq 1 then segmentplot,timdata,verdata,ye,last+1,2.,vcolor,6 segmentplot,timdata,hordata,ye,last+1,0,hcolor,6 endif oplot, [bt,et],[0.0,0.0], linestyle=2, color=1 ; ; Plot the keys to the horizontal and vertical velocity plot: plot, [0.0,1.0], [0.0,1.0], $ xstyle = 5, $ ystyle = 5, $ pos = [0.90,0.12,0.94,0.24], $ /nodata, color=1, /noerase ; Solid line for vertical: if vopt eq 1 then begin oplot, [0.0,1.0], [0.8,0.8], linestyle = 2, color = vcolor xyouts, 1.2, 0.75, 'VER', charsize=change, color=vcolor endif ; Solid line for horizontal (or horizontal minus corotation): oplot, [0.0,0.6], [0.4,0.4], linestyle = 0, color = hcolor if cor eq 0.0 then xyouts, 0.8, 0.35,'HOR-C', charsize=change,color=hcolor $ else xyouts, 1.4, 0.35,'HOR', charsize=change,color=hcolor RETURN END ;------------------------------------------------------------------------------ ; PRO segmnetplot - plot the line in segments so that there can be breaks ; in the line where time gaps in the data exist ;------------------------------------------------------------------------------ PRO segmentplot,tarray,yarray,ye,npts,lnstyle,icolor,samples_per_second ; ; tarray - array of x values ; yarray - array of y values ; ye - maximum range of yarray values that can be plotted = +/- ye ; npts - dimension of tarray and yarray ; lnstype - value for linestyple graphic key word ; icolor - index for color of line (exact color depends on palette) ; samples_per_second - self-explanatory ; xseg = fltarr(31) yseg = fltarr(31) deltat1 = 12./float(samples_per_second) ; Allow up to 12 samples to be missing ; before deciding that the line should be broken. For ; the driftmeter, this is 2 seconds of data. For SSM, ; this is 12 seconds of data (1 sec avg). ii = 0 NEXTSEG: if ii lt npts-1 then begin xseg(0) = tarray(ii) if yarray[ii] lt -9000. then yseg[0]=!VALUES.F_NAN else yseg[0]=yarray[ii] for jj=1,30 do begin if ii lt npts-1 then begin if tarray(ii+1)-tarray(ii) lt deltat1 then ii=ii+1 endif xseg(jj) = tarray(ii) if yarray[ii] lt -9000. then yseg[jj]=!VALUES.F_NAN else yseg[jj]=yarray[ii] endfor ;print,'seg_plot',xseg[0],yseg[0] oplot,xseg,yseg, linestyle=lnstyle, color=icolor, max_value=ye if ii lt npts-1 then begin if tarray(ii+1)-tarray(ii) ge deltat1 then ii=ii+1 endif if ii ge npts-1 then goto,HPLOTDONE goto,NEXTSEG endif HPLOTDONE: end ;----------------------------------------------------------------------------- ; ; procedure plot_ssm ; ; Sets up ranges for the magnetometer data, then plots it. This is ; a little more difficult than the other data, since it has to be plotted ; with meaningful times. Data gaps have been marked by setting the time values ; to -1000. This will cause problems if not handled properly. So time values ; are searched to find where there are negative and positive steps. ; Positive steps mean that there is really data, negative steps mean ; that the end of the data has arrived or there is a data gap, since ; time progresses only positively and -1000. implies data missing. ; So positive steps are plotted and negative steps are ignored. ; ; After the data is plotted, the ephemeris data is displayed. ; ;----------------------------------------------------------------------------- PRO plot_ssm,ssmdata, bt, et, nticks, tickval, $ yb, ye, ytick, ytickval, yticklab, $ begtime, nossm, ycol, $ change, output_type, dm_or_ssm_option ; ; Inputs; ; ssmdata- data array 3x1860(1 sample/sec = 31 min), ; 1) Time (sec), 2) Delta By, 3) Delta Bz ; bt - time of beginning of plot (sec) = 0. ; et - time of end of plot (sec) = 60*length of plot in minutes ; nticks - number of tick marks for X-axis ; tickval- values associated with tickmarks for X-axis ; yb - minimum values for Y-axis ; ye - maximum for Y-axis ; ytick - number of tick marks for Y-axis ; ytickval-values associated with tickmarks for Y-axis ; yticklab- labels to be placed at each Y-axis tick mark ; begtime - year,month,day,hour,minute,julian day (all floating point) ; for start of plot ; nossm - Flag to indicate whether SSM data file was opened (=0) or not (=1) ; ycol - Flag to indicate if the plot is to be done in color (=1) or B&W (=0) ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait ; output_type - index number for output device ; dm_or_ssm_option - Flag to indicate which data set are to be plotted ; =2, plot SSM data below J4/J5 spectrograms; ; =5 plot SSM data above J4/J5 color spectrograms with SM Ni data ; =4 plot SSM data alone above J4/J5 color spectrograms tickfake=strarr(20) for i=0,19 do tickfake(i)=' ' if dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin if dm_or_ssm_option eq 4 then begin ymin_ssm = 0.720 ymax_ssm = 0.820 endif else begin ymin_ssm = 0.720 ymax_ssm = 0.920 endelse endif else begin ymin_ssm = 0.120 ymax_ssm = 0.315 endelse ; Draw empty box for magnetometer data plot, [bt,et], [yb,ye], $ xstyle = 1, $ xticks = nticks, $ xtickname = tickfake, $ xtickv = tickval, $ xticklen = -0.03, $ ystyle = 1, $ yticks = ytick, $ ytickname = yticklab, $ ytickv = ytickval, $ yticklen = -0.005, $ yminor = 5, $ charsize = change, $ pos = [0.09,ymin_ssm,0.88,ymax_ssm], $ /noerase, color=1, $ /nodata if ycol eq 1.0 then begin ycolor = 40 ; color mode zcolor = 150 endif else begin ; B&W mode ycolor = 1 zcolor = 1 endelse if nossm eq 0 then begin ; magnetometer data file was opened, plot data utoff = begtime(3)*3600.0+begtime(4)*60.0 ; get timesteps: ssmsize = size(ssmdata) ssmd2 = ssmsize[2] dumdat = ssmdata[0,1:ssmd2-1] - ssmdata[0,0:ssmd2-2] ; where are they positve: nptsarr = where(dumdat gt 0.0,count) ; nptsarr=index of non-zero differences ;from one entry of ssmdata to the next print,size(nptsarr) ssmdatt = fltarr(count) ssmdaty = fltarr(count) ssmdatz = fltarr(count) for i=0,count-1 do begin k = nptsarr[i] ssmdatt[i] = ssmdata[0,k]-utoff ssmdaty[i] = ssmdata[1,k] ssmdatz[i] = ssmdata[2,k] endfor segmentplot,ssmdatt,ssmdaty,ye,count,1,ycolor,1 segmentplot,ssmdatt,ssmdatz,ye,count,0,zcolor,1 endif ; Put a dashed line down the y=0 line oplot, [bt,et],[0.0,0.0], linestyle=4,color=1 ; Plot the key to the y and z magnetometer component plot: plot, [0.0,1.0], [0.0,1.0], $ xstyle = 5, $ ystyle = 5, $ pos = [0.90,ymin_ssm+0.05*(ymax_ssm-ymin_ssm), 0.94, $ ymin_ssm+0.65*(ymax_ssm-ymin_ssm)], $ /nodata, color=1, $ /noerase ; Dotted line for y-data; oplot, [0.0,1.0], [0.8,0.8], linestyle = 4, color = ycolor !p.font=-1 xyouts, 1.3, 0.8, '!7D!6B!IY!N', charsize=change, color=1 ; Solid line for z-data; oplot, [0.0,1.0], [0.4,0.4], linestyle = 0, color = zcolor xyouts, 1.3, 0.4,'!7D!6B!IZ!N', charsize=change,color=1 if abs(output_type) lt 2. then !p.font=0 RETURN END ;----------------------------------------------------------------------------- ; ; procedure plot_smni - Plot the SSIES/SM data ; ;----------------------------------------------------------------------------- PRO PLOT_SMNI,smnidata, bt, et, nticks, tickval, ticklab, $ begtime, nosmni, change, ycol ; begtime - year,month,day,hour,minute,julian day (all floating point) ; for start of plot ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait ; nosmni - indicat whether SM/Ni data exists, 0=Yes; 1=No tickfake=strarr(20) for i=0,19 do tickfake(i)=' ' ymin_smni = 0.820 ymax_smni = 0.920 ; find maximum value of Log10[ Ni ] if nosmni eq 0 then maxlogni = max(smnidata[1,*]) else maxlogni = 5.5 print,'Max Log10 Ni =', maxlogni if maxlogni gt 5.0 then begin ymin_value = 4.0 ymax_value = 6.0 endif else begin ymin_value = 3.0 ymax_value = 5.0 endelse plot, [bt,et], [10.^ymin_value,10.^ymax_value], $ xstyle = 1, $ xticks = nticks, $ xtickname = tickfake, $ xtickv = tickval, $ xticklen = -0.03, $ /ylog, ystyle = 1, $ yticklen = -0.008, $ ytickname = tickfake, $ charsize = change, $ pos = [0.09,ymin_smni,0.88,ymax_smni], $ /noerase, color=1, $ /nodata if ycol eq 1.0 then begin ycolor = 80 ; color mode zcolor = 150 endif else begin ; B&W mode ycolor = 1 zcolor = 1 endelse for i=0,2 do begin st = string(format='(a4,i1,a2)','10!E',ymin_value+i,'!N') xyouts,et+0.02*(et-bt),10^(ymin_value+i),st,color=1 endfor utoff = begtime(3)*3600.0+begtime(4)*60.0 if nosmni ne 0 then goto,DONE_PLOT_SMNI ; get timesteps: smnisize = size(smnidata) smnid2 = smnisize[2] dumdat = smnidata(0,1:smnid2-1) - smnidata(0,0:smnid2-2) ; where are they positve: nptsarr = where(dumdat gt 0.0,count) if count lt 1 then goto,DONE_PLOT_SMNI dumdat = nptsarr(1:count-1) - nptsarr(0:count-2) ; If there are data gaps greater than 1 second, where are they: poss = where(dumdat gt 1.0,count) ; COUNT is the number of good plot ranges - 1, so if there are ; no data gaps, COUNT = 0 print,'Count= ',count if count gt 0 then begin for i=0,count do begin ; If it's the first time range, the first data point will be the ; first element in the array, if it's not, the first element will ; be the first good point after the last good point in the previous ; good data range: if i eq 0 then first = 0 else first = nptsarr(poss(i-1)+1)+1 ; If it's the last time range, the last point will be the last good ; point, if not, it will be last good element in the current time ; range: if i lt count then last = nptsarr(poss(i))+1 $ else last = nptsarr(n_elements(nptsarr)-1)+1 ; Grab the good time range: timdata = smnidata(0,first:last)-utoff ydata = smnidata(1,first:last) ii2 = last-first for ii=0,ii2 do begin if ydata[ii] le 0. then ydata[ii] =!VALUES.F_NAN else $ ydata[ii]=10.^ydata[ii] ; Convert from Log10 Ni to Ni endfor timdum = timdata-float(fix(timdata)) srch = where(timdum eq 0.0, scount) ; Plot good data: oplot, timdata, ydata, max_value=10.^ymax_value, linestyle=0, color=ycolor endfor endif else begin ; If all the data is good, decide where the end of the data is, ; grab it and plot it: last = nptsarr(n_elements(nptsarr)-1)+1 timdata = smnidata(0,0:last)-utoff ydata = smnidata(1,0:last) for i=0,last do begin ; Convert from Log10 Ni to Ni if ydata[i] gt 0. then ydata[i]=10.^ydata[i] else ydata[i]=!VALUES.F_NAN endfor timdum = timdata-float(fix(timdata)) srch = where(timdum eq 0.0, scount) oplot, timdata, ydata, max_value=ye, linestyle=0, color=ycolor endelse DONE_PLOT_SMNI: RETURN END ;----------------------------------------------------------------------------- ; ; procedure plot_ee ; ; This procedure plots out the total energy flux and the average energy ; flux for the time period selected. The program filters the data if ; that is what the user wants (calls RAVE). ; ;----------------------------------------------------------------------------- PRO plot_ee, eeflx, eave, ieflx, iave, ticklab, intsw, $ nticks, tickval, bt, et, yfil, yticklab, change ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait ; Blank out every other tick name with a space ('' does not cause any ; change to the tick name) ytickn = [' ',' ',' ','',' ','',' ','',' ',' ',' ',' '] if (intsw eq 1) then begin ; Set up plotting range for the energy flux. This range was used, since ; it includes all of the ranges specified by the Newell and Meng papers ; on identification of mapped regions. plot, [bt,et], [10.0^9.0,10.0^13.0], $ /ylog, xstyle = 1, $ xticks = nticks, $ xtickname = ticklab, $ xtickv = tickval, $ xticklen = 0.03, $ yminor = 0, $ ytickname = ytickn, $ yticklen = -0.005, $ charsize = change, $ pos = [0.09,0.820,0.88,0.920], $ /noerase, color=1, $ /nodata ; Filter the data if the user wants: n_filter = 9 if yfil eq 0.0 then rave, eeflx, n_filter, stdev if yfil eq 0.0 then rave, ieflx, n_filter, stdev ; Plot the data: oplot, eeflx, max_value = 10.0^20.0, color = 1 oplot, ieflx, max_value = 10.0^20.0, linestyle = 1, color = 1 endif if (intsw eq 0) then begin ; Set up plotting range for the number flux. This range was used, since ; it includes all of the ranges specified by the Newell and Meng papers ; on identification of mapped regions. ; plot, [bt,et], [10.0^5.0,10.0^9.0], $ /ylog, xstyle = 1, $ xticks = nticks, $ xtickname = ticklab, $ xtickv = tickval, $ xticklen = 0.03, $ yminor = 0, $ ytickname = ytickn, $ yticklen = -0.005, $ charsize = change, $ pos = [0.09,0.820,0.88,0.920], $ /noerase, color=1, $ /nodata ; Filter the data if the user wants: n_filter = 9 if yfil eq 0.0 then rave, eeflx, n_filter, stdev if yfil eq 0.0 then rave, ieflx, n_filter, stdev ; Plot the data: oplot, eeflx, max_value = 10.0^16.0, color=1 oplot, ieflx, max_value = 10.0^16.0, color=1, linestyle = 1 endif ; Set up new y tick names for the average energy: ytickn = ['10!A2','10!A3','10!A4'] ; Set up the plotting range (same as the image ranges): nt1 = nticks + 1 plot, [bt,et], [31.0,30500.0], $ /ylog, xstyle = 1, $ ystyle = 1, $ xticks = nticks, $ xtickname = ticklab, $ xtickv = tickval, $ xticklen = 0.03, $ yminor = 0, $ ytickname = ytickn, $ yticklen = -0.005, $ charsize = change, $ pos = [0.09,.720,0.88,0.820], $ /noerase, color=1, $ /nodata ; Filter the data if the user wants: n_filter = 9 if yfil eq 0.0 then rave, eave, n_filter, stdev if yfil eq 0.0 then rave, iave, n_filter, stdev ; Plot the data: oplot, eave, max_value = 10.0^20.0, color=1 oplot, iave, max_value = 10.0^20.0, color=1, linestyle = 1.0 ; Set up the key to the plot: plot, [0.0,1.0], [0.0,1.0], $ xstyle = 5, $ ystyle = 5, $ pos = [0.91,0.70,0.94,0.90], $ /nodata, color=1, $ /noerase ; Ions with dotted lines: oplot, [0.0,1.0], [0.7,0.7], linestyle = 1 xyouts, 1.5, 0.7, 'Ion', charsize=change,color=1 ; electrons with solid lines: oplot, [0.0,1.0], [0.5,0.5], linestyle = 0 xyouts, 1.5, 0.5, 'Ele', charsize=change,color=1 RETURN END ;----------------------------------------------------------------------------- ; ; procedure plot_color_bar ; ; plots the color (or black and white) table for the particle color spectrogram ; on the side of the plot ; ;----------------------------------------------------------------------------- PRO plot_color_bar, cyb, cye, ctick, ctickval, cticke, cticki, change, $ diswitch, output_type, dm_or_ssm_option,ycol ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait ; diswitch - flag to indicate whether eflx and iflx are differential ; number flux (=0) or differential energy flux (=1) ; dm_or_ssm_option - flag to indicate whether DM data (=1), SSM data (=2), ; neither (=3) or both (=5) or both plus SM Ni data (=5) ; will be plotted with the particle data ; YCOL - flag to indicate color (= 1.0) or b/w (= 0.0) plot ; cticke = ['',' ','',' ','',' '] ; Set up plot range and position. We don't want ticks on the bottom ; of the color chart, since there is no x values, so we plot with a ; xstyle=5 which is not x axis. Afterwards, we put in a blank axis. if diswitch eq 0 then begin ;electron diff. number flux top = 7.0 bot = 1.68 endif else begin ;electron diff. energy flux top = 10.0 bot = 4.68 endelse if dm_or_ssm_option ne 3 then begin ymin_bar = 0.521 ymax_bar = 0.720 endif else begin ymin_bar = 0.421 ymax_bar = 0.720 endelse ; ; Preset the loop value to 200 lpval = 200 ; plot, [0.0,1.0],[10.0^(bot),10.0^(top)], $ /ylog, xstyle = 5, $ ystyle = 1, $ yminor = 0, $ ytickname = cticke, $ yticklen = -0.2, $ pos = [0.93,ymin_bar,0.95,ymax_bar], $ charsize = change, $ /nodata, color=1, $ /noerase ; The blank axis (top and bottom): oplot, [0.0,1.0],[10.0^(bot),10.0^(bot)] oplot, [0.0,1.0],[10.0^(top),10.0^(top)] ; Next, we set up a plot range, just inside the plotted axis, so when we ; plot the colors, the black outline will not be erased: plot, [0.0,1.0],[cyb,cye], $ xstyle = 5, $ ystyle = 5, $ pos = [0.932,ymin_bar + 0.001,0.949,ymax_bar - 0.001], $ charsize = change, $ /nodata, color=1, $ /noerase ; Plot the 200 (or 60) colors using POLYFILLs. This garantees that no matter how ; big or small the color table is, it will always have 200 distinct, ; touching color values. This is much better than just OPLOTting with ; different line thickness statements. for i = 1, lpval do begin polyx = [0.0,0.0,1.0,1.0,0.0] sy = cyb+float(i-1)*(cye-cyb)/float(lpval) ey = cyb+float(i)*(cye-cyb)/float(lpval) polyy = [sy,ey,ey,sy,sy] if ycol gt 0 then polyfill, polyx, polyy, color = i else $ polyfill, polyx, polyy, color = i+40 endfor ; Label the color table: xyouts, 2.0, cyb+0.5*(cye-cyb),'ELECTRON FLUX', color=1, $ alignment=0.5, orientation = 90.0, charsize=change*0.9 if diswitch eq 0 then begin xyouts, 3.0, cyb+0.5*(cye-cyb),'particles/(cm!E2!N s sr eV)', color=1, $ alignment=0.5, orientation = 90.0, charsize=change*0.9 endif else begin xyouts, 3.0, cyb+0.5*(cye-cyb),'eV/(cm!E2!N s sr eV)', color=1, $ alignment=0.5, orientation = 90.0, charsize=change*0.9 endelse ; Do the exact same thing with the ion color table, but with a different ; range of values: cticke = ['',' ','',' ','',' '] if diswitch eq 0 then begin top = 5.0 bot = -0.32 endif else begin top = 8.0 bot = 2.68 endelse if dm_or_ssm_option ne 3 then begin ymin_bar = 0.320 ymax_bar = 0.510 endif else begin ymin_bar = 0.120 ymax_bar = 0.410 endelse plot, [0.0,1.0],[10.0^(bot),10.0^(top)], $ /ylog, xstyle = 5, $ ystyle = 1, $ yminor = 0, $ ytickname = cticke, $ yticklen = -0.2, $ pos = [0.930,ymin_bar,0.950,ymax_bar], $ charsize = change, $ /nodata, color=1, $ /noerase oplot, [0.0,1.0],[10.0^(bot),10.0^(bot)] oplot, [0.0,1.0],[10.0^(top),10.0^(top)] plot, [0.0,1.0],[cyb,cye], $ xstyle = 5, $ ystyle = 5, $ yticklen = -0.1, $ pos = [0.932,ymin_bar+0.001,0.949,ymax_bar-0.001], $ charsize = change, $ /nodata, color=1, $ /noerase for i = 1, lpval do begin polyx = [0.0,0.0,1.0,1.0,0.0] sy = cyb+float(i-1)*(cye-cyb)/float(lpval) ey = cyb+float(i)*(cye-cyb)/float(lpval) polyy = [sy,ey,ey,sy,sy] if ycol gt 0 then polyfill, polyx, polyy, color = i else $ polyfill, polyx, polyy, color = i+40 endfor xyouts, 2.0, cyb+0.5*(cye-cyb),'ION FLUX', color=1, $ alignment=0.5, orientation = 90.0, charsize=change*0.9 if diswitch eq 0 then begin xyouts, 3.0, cyb+0.5*(cye-cyb),'particles/(cm!E2!N s sr eV)', color=1, $ alignment=0.5, orientation = 90.0, charsize=change*0.9 endif else begin xyouts, 3.0, cyb+0.5*(cye-cyb),'eV/(cm!E2!N s sr eV)', color=1, $ alignment=0.5, orientation = 90.0, charsize=change*0.9 endelse RETURN END ;----------------------------------------------------------------------------- ; ; procedure plot_extra ; ; This procedure simply plots a couple of extra things on the page, such ; as the satellite number, the date,and the plot labels ; ;----------------------------------------------------------------------------- PRO plot_extra, begtime, monname, satel, change, psfile, dm_or_ssm_option, $ bt, et, nticks, tickval, ticklab, imul, mlt, mlat, $ yb, ye, intsw, output_type, removebackground, JGF_LABEL ; begtime - year,month,day,hour,minute,julian day (all floating point) ; for start of plot ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait tickfake = strarr(20) dayofmon = strcompress(string(begtime(2)),/remove_all) year = strcompress(string(begtime(0)),/remove_all) month = strmid(monname,3*(begtime(1)-1),3) ; date = dayofmon + ' ' + month + ' ' + year fsat = 'F' + strcompress(string(satel),/remove_all) if dm_or_ssm_option eq 1 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then $ unit1 = 'M/SEC' if dm_or_ssm_option eq 2 then unit1 = ' nT ' if dm_or_ssm_option eq 3 then unit1 = ' ' unit2 = '(eV)' unit3 = 'Electrons' unit4 = 'Ions' if dm_or_ssm_option ne 4 and dm_or_ssm_option ne 5 then begin if (intsw eq 1) then unit5 = '!7U!6!IE!N!X' if (intsw eq 0) then unit5 = '!7U!6!IN!N!X' unit6 = '!6E!Iave!N' endif else begin unit5 = 'N!Ii!N (cm!E-3!N)' unit6 = 'nT' endelse if (intsw eq 1) then unit7 = 'eV/(cm!E2!N s sr)' if (intsw eq 0) then unit7 = ' 1/(cm!E2!N s sr)' unit8 = 'eV' xyouts, 0.02, 0.95, fsat, charsize = 2.0*change,color=1, /norm xyouts, 0.04, 0.22, unit1, alignment = 0.5, orientation = 90.0,color=1, $ charsize=change, /norm xyouts, 0.5, 0.95, date, alignment = 0.5, charsize = 2.0*change,color=1, /norm if removebackground NE 0 then xyouts,0.5,0.93,'SSJ Background Removed',$ alignment=0.5, charsize = 0.75*change,color=1, /norm if dm_or_ssm_option eq 3 then ypos =0.42 else ypos = 0.52 xyouts, 0.04, ypos, unit2, alignment = 0.5, orientation = 90.0, $ charsize=change,color=1, /norm xyouts, 0.01, ypos, JGF_LABEL, alignment = 0.5, orientation = 90.0, $ charsize=change*0.5,color=1, /norm print,'plot_extra: JGF_LABEL = ', JGF_LABEL if dm_or_ssm_option eq 3 then ypos =0.570 else ypos = 0.620 xyouts, 0.04, ypos, unit3, alignment = 0.5, orientation = 90.0, $ charsize=change,color=1, /norm if dm_or_ssm_option eq 3 then ypos =0.27 else ypos = 0.42 xyouts, 0.04, ypos, unit4, alignment = 0.5, orientation = 90.0, $ charsize=change,color=1, /norm if dm_or_ssm_option ne 4 and dm_or_ssm_option ne 5 then begin !p.font = -1 ; vector drawn fonts to be used xyouts, 0.02, 0.870, unit5, /norm, alignment=0.5, orientation = 90.0, $ color=1,charsize = 1.5*change xyouts, 0.02, 0.773, unit6, /norm, alignment=0.5, orientation = 90.0, $ color=1,charsize = 1.5*change endif else begin if dm_or_ssm_option eq 4 then begin xyouts, 0.04, 0.870, unit5, /norm, alignment=0.5, orientation = 90.0, $ color=1,charsize = change xyouts, 0.04, 0.773, unit6, /norm, alignment=0.5, orientation = 90.0, $ color=1,charsize = change endif else begin xyouts, 0.04, 0.820, unit6, /norm, alignment=0.5, orientation = 90.0, $ color=1,charsize = change endelse endelse if abs(output_type) le 2 then !p.font = 0 ; Graphics device specific font to be used ; /times only works for a postscript file, not for an X window if psfile ne '' then begin ; set for Time-Roman postscript font device, /times endif if dm_or_ssm_option ne 4 and dm_or_ssm_option ne 5 then begin xyouts, 0.05, 0.770, unit8, /norm, alignment=0.5, orientation = 90.0, $ color=1,charsize = change*0.8 xyouts, 0.05, 0.870, unit7, /norm, alignment=0.5, orientation = 90.0, $ color=1,charsize = change*0.8 endif ; ; This defines the area for the ephemeris plot, [bt,et], [yb,ye], $ xstyle = 5, $ ystyle = 5, $ pos = [0.09,0.120,0.88,0.315], $ /noerase, color=1, $ /nodata ; ; Plot Ephemeris data. In order to do this so the data ends up in the ; right place every time (vertically), a normalized (to range) ; coordinate system is used. ; This is done with a series of XYOUT commands to output the text to the ; specific tick mark locations and offset vertically for each ; different piece of data (such as lat, mlt) ; range = ye - yb for i=0,nticks do begin utmin = fix(begtime(4) + tickval(i)/60.0) uthr = begtime(3) if utmin ge 60 then begin utmin = utmin - 60 uthr = uthr + 1 endif if utmin ge 10 then $ utstr = strcompress(string(uthr)+':'+string(utmin),/remove_all) $ else $ utstr = strcompress(string(uthr)+':0'+string(utmin),/remove_all) ; Since we have the lat and mlt at each minute, but only need it at ; specific tick marks, IMUL*I is used as a counter to pick out each ; tick mark time.... so if we had 18 minutes worth of data and 9 tick ; marks, IMUL would be 2... mlatstr = strcompress(string(mlat(i*imul)),/remove_all) if mlat(i*imul) gt 0.0 then mlatstr = strmid(mlatstr,0,4) $ else mlatstr = strmid(mlatstr,0,5) mltstr = strcompress(string(mlt(i*imul)),/remove_all) mltstr = strmid(mltstr,0,4) ; Since we use range as a coordinate system, the position of the tick ; marks will be the same every time (vertically). yl = yb-0.25*range xyouts, tickval(i), yb-0.14*range, ticklab(i), alignment=0.5, $ color=1,charsize=change xyouts, tickval(i), yb-0.28*range, utstr, alignment=0.5, $ color=1,charsize=change if abs(mlat(i*imul)) eq 0.0 and abs(mlt(i*imul)) eq 0.0 then GOTO,SKIP_EPH xyouts, tickval(i), yb-0.42*range, mlatstr, alignment=0.5, $ color=1,charsize=change xyouts, tickval(i), yb-0.56*range, mltstr, alignment=0.5, $ color=1,charsize=change SKIP_EPH: endfor ; put the labels for the ephemeris data under the color ; table, but also below the DM data. plot, [et,et+1.0],[yb,ye], $ xstyle = 5, $ ystyle = 5, $ pos = [0.96,0.12,0.98,0.32], $ /nodata, color=1, $ /noerase xyouts, et, yb-0.14*range,' UT(SEC)', alignment=0.5,color=1, charsize=change xyouts, et, yb-0.28*range,' HH:MM', alignment=0.5,color=1, charsize=change xyouts, et, yb-0.42*range,' MLAT', alignment=0.5,color=1, charsize=change xyouts, et, yb-0.56*range,' MLT', alignment=0.5,color=1, charsize=change RETURN END ;----------------------------------------------------------------------------- ; ; procedure add_points ; ; This procedure allows the user to add vertical lines, labels for the ; vertical lines and points to the output. The plot does not look ; very nice when it is output to the display, but the resulting ps ; file will look exactly the same as usual and with extra lines and ; points. ; The reasonning for this type of procedure is for the identification and ; labeling of mapped regions or oddities in the data. ; PRO add_points, vlines, pts, bt, et, vcount, pcount,$ vstr, vlabx, vlabc, change, baset, datafile_or_keyboard_input, $ dm_or_ssm_option ; vcount - Number of vertical lines ; vlines - X value for vertical lines ; vlabc - Number of labels ; vstr - String for each label. Array size expands as user enters labels. ; vlabx - X value for each label ; pcount - Number of symbols to be added to filled circle symbol(s) ; pts - X,Y coordinate of position(s) to place filled cirle symbol(s) ; (Dec 2002. The program does not have any code to get input from ; the user for placement of these symbols. Until the code is added ; this option is non-functional.) ; bt - beginning time of plot (always zeo) ; et - end time of plot (=number of minute in plot times 60) ; baset - time of beginning of plot in seconds of the day ; change - scaling factor for lettering and characters on plot. Nominally ; 1.0 for landscape and 0.75 for portrait ; datafile_or_keyboard_input ; if 0, read input from keyboard ; if 1, readf input from a prepared file ; dm_or_ssm_option - if 1, plot dm data with SSJ4 or SSJ5 data ; if 2, plot SSM data with SSJ4 or SSJ5 data; ; if 3, plot only SSJ4 or SSJ5 data ; if 4, plot DM, SSM and SM Ni data with SSJ4/5 spectrogram ; if 5, plot DM, and SSM data with SSJ4/5 spectrogram if n_elements(pcount) eq 0 then pcount = 0 if n_elements(pts) eq 0 then pts = fltarr(1) ; Create a filled circle symbol to use as a marker point (this filled ; circle is NOT a basic IDL symbol, that is why we have to create it). ; A user defined symbol is used by plotting with a psym=8. a = findgen(30)*2.0*!pi/29.0 usersym, 0.5*sin(a),0.5*cos(a),/fill ; set up the plotting coordinates, but over the data which already ; exists plot, [bt,et],[0.0,1.0], pos = [0.09,0.12,0.88,0.92], $ xstyle=5, ystyle=5, /nodata, /noerase,color = 1 ; These three statements have to be made even if the user is using a ; postscript file, since they could have added points, then selected to ; get a hard copy of the plot. If this occurs, the plotting range and ; symbols need to be set up so the hard copy can have the lines and ; points. ; Only ask about points and lines if the user has selected to NOT plot to ; a postscript file ; ; if datafile_or_keyboard_input eq 0 and dm_or_ssm_option eq 1 then begin ; endif else begin if vcount gt 0 then $ for i=0,vcount-1 do begin oplot,[vlines[i]-baset,vlines[i]-baset],[0.0,1.0], thick=2.0 endfor if pcount gt 0 then $ for i=0,pcount-1 do $ oplot, [pts(i,0)],[pts(i,1)], psym=8 if vlabc gt 0 then $ for i=0,vlabc-1 do $ xyouts, vlabx(i)-baset, 1.005, vstr(i), $ alignment = 0.5,color=1, charsize=0.7*change ; endelse return end ;----------------------------------------------------------------------------- ; ; procedure rave ; ; this procedure simply averages 9 points together and then calls the midpoint ; the average. So is is basically a low pass filter, doing a running average. ; I need to pass another variable to the procedure to make the 9 points a ; general number of points. ; ;----------------------------------------------------------------------------- PRO rave, avearr, nadd, stdev ; How many points are there in the array: npts = n_elements(avearr)-1 stdev = fltarr(npts+1) ; If there is enough data, go ahead with the filter if (npts gt nadd) and (nadd ge 3) then begin ; We need another array, so we don't double count points. In other words, ; filter the old array, while shoving the filtered values into a new ; array. dumarr = fltarr(npts+1) ; Since the first 4 and last 4 points are not filtered, just put them ; into the new array left = (nadd-1)/2 right = nadd - left - 1 dumarr(0:left-1) = avearr(0:left-1) dumarr(npts-right-1:npts) = avearr(npts-right-1:npts) ; Filter the rest for i = left, npts-right do begin sum = 0.0 for j=-left,right do sum=sum+avearr(i+j) dumarr(i) = sum / nadd sum = 0.0 for j=-left,right do sum=sum+abs(dumarr(i)-avearr(i+j)) stdev(i) = sum / nadd endfor ; Put the new data back into the original data array avearr = dumarr endif RETURN END ;----------------------------------------------------------------------------- ; ; procedure interpo ; ; This procedure creates the bitmap that is used to show the J4/J5 electron/ion ; spectograms. Originally, this bitmap was only valid Postscript and HP plot ; becuase the device "stretched" the bitmap to fit the allow space. For the ; terminal modes, the bitmap must exactly fit the allowed space. ; ; For PS or HP plot, the x-dimension of the bitmap = number of seconds in plot ; and the y-dimension is created by changing the 20 energy bins up to 96 bins. ; For X or WIN plot, the x and y dimension are set by the window, The ; arrays of n_seconds by 20 energies must be expanded/shrunk to fit. ; After the map of fluxes have been created, then each pixel is assinged ; a "color" index (a number between 1 and 200). These colors are ; then corrected for bad values, such as negative numbers are assigned to ; be white, larger than the maximum color are assigned to the maximum, ect... ; ;----------------------------------------------------------------------------- PRO interpo, eflx, iflx, npts, ele2d, ion2d, diswitch, output_type, dm_or_ssm_option ; Inputs: ; eflx - array of 60*npts by 20 values of electron number or energy flux ; iflx - array of 60*npts by 20 values of ion number or energy flux ; npts - number of minutes in plot. ; diswitch - flag to indicate whether eflx and iflx are differential ; number flux (=0) or differential energy flux (=1) ; output_type - flag to indicate destination device for plot ; Outputs: ; ele2d - bitmap of electron spectogram ; ion2d - bitmap of ion spectogram ; Set the base value (10^1.68 for electron number, ; 10^4.68 for electron energy). We want to subtract off the right value. if diswitch eq 0 then begin ebase = 1.68 ibase = -0.32 endif else begin ebase = 4.68 ibase = 2.68 endelse ;-------------------------------------------------- ; If PS or HP plot then make and file fixed arrays if abs(output_type) le 2 then begin ; Initialize arrays ele2d = fltarr(npts*60,91) ion2d = fltarr(npts*60,91) ; Loop through old array, while filling the new, larger, array with ; interpolated values. The old array is backwards, so the new array ; looks like it is being filled from the top down (which it is). for i=0,8 do for j=0,4 do begin ele2d(*,90-(i*5+j)) = $ float(j)*(eflx(*,i+1)-eflx(*,i))/5.0 + eflx(*,i) ion2d(*,90-(i*5+j)) = $ float(j)*(iflx(*,i+1)-iflx(*,i))/5.0 + iflx(*,i) endfor ; Skip Channel 11. For SSJ4, Channel 11 is redundant with 10. For ; SSJ5, Channel 11 does not contain particle data. i=9 for j=0,4 do begin ele2d(*,90-(i*5+j)) = $ float(j)*(eflx(*,i+2)-eflx(*,i))/5.0 + eflx(*,i) ion2d(*,90-(i*5+j)) = $ float(j)*(iflx(*,i+2)-iflx(*,i))/5.0 + iflx(*,i) endfor for i=11,18 do for j=0,4 do begin ele2d(*,90-(i*5+j-5)) = $ float(j)*(eflx(*,i+1)-eflx(*,i))/5.0 + eflx(*,i) ion2d(*,90-(i*5+j-5)) = $ float(j)*(iflx(*,i+1)-iflx(*,i))/5.0 + iflx(*,i) endfor ; The last value is just stuck in ele2d[*,0] = eflx[*,19] ion2d[*,0] = iflx[*,19] ; End of PS or HP plot ( abs(output_type) = 1 or 2 ) endif else begin ; Start creation of arrays for terminal plots ( abs(output_type) = 3 ) ; The Window size of 750 by 650 was set by the Window command in ; the MAIN procedure just after the set_plot command. nx = fix(750*0.789) if dm_or_ssm_option ne 3 then ny = fix(650*0.194) else ny = fix(650*0.294) ; Initialize arrays ele2d = fltarr(nx,ny) ion2d = fltarr(nx,ny) tempe = fltarr(20) tempi = fltarr(20) ele2d[*,*] = 1000. ion2d[*,*] = 100. ; First expand/contract the input array in the x-direction to fit the ; window. deltax = float(60*npts)/nx xstart = 0. ixstart = fix(xstart) xend = deltax ixend = fix(xend) ixnew = 0 while ixend lt 60*npts do begin for i=0,18 do begin ii = 18-i ik = i if ik ge 10 then ik = ik + 1 ; Skip Channel 11. It is redundatn in SSJ4 and not ; meaningful in SSJ5 sume = 0. sumi = 0. nume = 0 numi = 0 for j=ixstart,ixend do begin if eflx(j,ik) gt 0 then begin sume = sume + eflx(j,ik) nume = nume + 1 endif if iflx(j,ik) gt 0 then begin sumi = sumi + iflx(j,ik) numi = numi + 1 endif endfor if nume gt 0 then ele2d(ixnew,ii) = sume/float(nume) else ele2d(ixnew,ii) = 0. if numi gt 0 then ion2d(ixnew,ii) = sumi/float(numi) else ion2d(ixnew,ii) = 0. endfor xstart = xstart + deltax ixstart = fix(xstart) xend = xend + deltax ixend = fix(xend) ixnew = ixnew + 1 endwhile ; Now expand the array in the y direction deltay = float(ny-1)/18. for j=0,nx-1 do begin ystart = 0. iystart = fix(ystart) yend = deltay iyend = fix(yend) for i=0,18 do begin tempe(i) = ele2d(j,i) tempi(i) = ion2d(j,i) endfor for i=0,17 do begin for ii=iystart,iyend do begin ele2d(j,ii) = tempe(i)+ (tempe(i+1)-tempe(i))*float(ii-iystart)/deltay ion2d(j,ii) = tempi(i)+ (tempi(i+1)-tempi(i))*float(ii-iystart)/deltay endfor ystart = ystart + deltay iystart = fix(ystart) yend = yend + deltay iyend = fix(yend) endfor endfor ; End of X or MS Window Plot ( abs(output_type) = 3 ) endelse ; Assign a color for each electron value which is greater than 0 ; electron base is ebase (see above) logarr = where(ele2d gt 0.0, count) if count gt 0 then $ ele2d(logarr) = 40.0*(alog(ele2d(logarr))/2.30259 - ebase) ; Check for "errors" logarr = where(ele2d gt 200.0, count) if count gt 0 then ele2d(logarr) = 200.0 logarr = where(ele2d lt 1.0, count) if count gt 0 then begin ; no counts = background color = white if abs(output_type) le 2 then ele2d(logarr) = 255. else ele2d(logarr) = 0. endif ; Do the same for the ions, but with a different base (ibase, see above) logarr = where(ion2d gt 0.0, count) if count gt 0 then $ ion2d(logarr) = 40.0*(alog(ion2d(logarr))/2.30259 - ibase) logarr = where(ion2d gt 200.0, count) if count gt 0 then ion2d(logarr) = 200.0 logarr = where(ion2d lt 1.0, count) if count gt 0 then begin ; no counts = background color = white if abs(output_type) le 2 then ion2d(logarr) = 255. else ion2d(logarr) = 0. endif RETURN END ;----------------------------------------------------------------------------- ; ; procedure calc_flx ; ; This section gets the conversion factors from electron and ion ; counts to number flux and energy flux for the different energy gates. ; and computes the integrated fluxes, differential fluxes, avg energy ; and fluxes from 450 eV to 30 keV. ; ; ; Variables Used - ; ecnt, icnt - counts from J4 or J5 instrument. ; eflx, iflx - differential directional number flux in units of ; particles/(cm^2 s sr eV). ; eeflx, ieflx - integral energy flux in units of eV/(cm^2 s sr) ; printed out in units of mW/m^2 (or ergs/cm^2-s) with fac=1.602e-12*pi ; 1eV = 1.602e-12 erg and pi/sr (ster-radian) to get the isotropic downgoing particles ; 2piR is circumference of circle, so half circumference is piR. Thus pi for downgoing only. ; mW/m^2 = 1.e-3J/cm^2*1.e+4-s = 1.e-7*1.e+7erg/cm^2-s = erg/cm^2-s ; NOTE: !pi is a special term for pi in IDL (! is NOT a comment or inverse etc) ; enflx, inflx - integral number flux in units of part/(cm^2 s sr) ; eave, iave - average energy in eV = eeflx/enflx or ieflx/inflx ; oeeflx, oieflx- energy flux into ionosphere assuming isotropic ; flux; for this, energies below ~450 eV are zeroed; ; units of eV/(cm^2 s). The energies which are zeroed ; are the 14-20 channels. ; oeave, oiave - average energy in eV for energies above and ; including ~450 eV. ; deeflx, deiflx- differential directional energy flux in units of ; eV/(cm^2 s sr eV). Calculated by multiplying eflx ; and iflx by the central energy of each bin and by ; the energy independent geometric factor, GF. The GF ; multiplication cancels the previous division by GF ; used to derive eflx and iflx. ; begtime - year,month,day,hour,minute,julian day (all floating point) ; for start of plot ;----------------------------------------------------------------------------- PRO calc_flx, ecnt, icnt, isat, eeflx, ieflx, eave, iave, eflx, iflx, $ oeeflx, oieflx, enflx, inflx, deeflx, deiflx, output_type, iyear, $ ; begtime, JGF_LABEL ; Emery 05/07: added begtime to calc_flx ; Emery 10/08: added minimum_cglat,mlat1m,mlt1m to calc_flx begtime,JGF_LABEL,totlen,minimum_cglat,uts1m,glat1m,glon1m,alt1m, $ glat1101m,glon1101m,mlon1101m,mlat1m,mlt1m ; Emery 10/08: added totlen, dumlat and ephemeris data to calc_flx (uts1m use float(hour)*3600., or bad!) ELECAL0=FLTARR(20) ; e energy of channel i ELECAL1=FLTARR(20) ; e CEFLX - Conversion factor to electron integral energy flux ELECAL2=FLTARR(20) ; e CNFLX - Conversion factor to electron integral number flux ELECAL3=FLTARR(20) ; e CFLX - Conversion factor to electron diff number flux ELEGF =FLTARR(20) ; e geometric factors IONCAL0=FLTARR(20) ; i energy of channel i IONCAL1=FLTARR(20) ; i CEFLX - Conversion factor to ion integral energy flux IONCAL2=FLTARR(20) ; i CNFLX - Conversion factor to ion integral number flux IONCAL3=FLTARR(20) ; i CFLX - Conversion factor to ion diff number flux IONGF =FLTARR(20) ; i geometric factors JGF_LABEL = ' ' ; text string indicating the date that the SSJ calibration was released read_geomfac, isat, ELECAL0, ELECAL1, ELECAL2, ELECAL3, $ ELEGF, IONCAL0, IONCAL1, IONCAL2, IONCAL3, $ IONGF, iyear, JGF_LABEL print,'calc_flx: JGF_LABEL = ', JGF_LABEL ; Create the arrarys for integrated number and energy fluxes fac = 1.602e-12*!pi eeflx450 = fltarr(n_elements(eflx(*,0))) enflx450 = fltarr(n_elements(eflx(*,0))) ieflx450 = fltarr(n_elements(eflx(*,0))) inflx450 = fltarr(n_elements(eflx(*,0))) eeflx = fltarr(n_elements(eflx(*,0))) ieflx = fltarr(n_elements(eflx(*,0))) enflx = fltarr(n_elements(eflx(*,0))) inflx = fltarr(n_elements(eflx(*,0))) ; This next section actually calculates the electron and ion integrated energy ; and number flux, using the calibration numbers just read in. ; Open file for printing integrated fluxes ; Kilcommons - This is the location in the code in which the ascii output is created. It only contains SSJ flux data, obviously. get_lun,u_eflux s = string(format='(a18,i2.2,a1,i4.4,4i2.2,a4)','cenergy_int_flux_f',isat,'_',begtime[0],begtime[1],begtime[2],begtime[3],begtime[4],'.txt') openw,u_eflux,s s = string(format='(a1,i2.2,a1,i4.4,4i2.2)','F',isat,'_',begtime[0],begtime[1],begtime[2],begtime[3],begtime[4]) ; format for 15sec DMSP data from Fred Rich read in sef15r.f to get AMIE input ; read (11,"(i5,i4,1x,i2,1x,i2,1x,i2,2x,3f7.1,1x,3f7.1,2x,2e11.3, ; | f7.2,1x,2e11.3,f7.2,i6)",end=500) iyr2,ida2,ihr2,imn2,isc2,glat2, ; | glon2,alt2,cmlat2,cmlon2,rmlt2,totjel2,totejel2,ekev2,totji2, ; | toteji2,ekevi2,ksec2 ;printf,u_eflux,s,' Pt. No, Int_ele_eflx, Int_ion_eflx, Int_ele_nflx, Int_ion_nflx' printf,u_eflux,s,' #,utsec,glat,glon,alt,glat110,glon110,mlat110,mlon110,mlt,Int_ele_eflx(mW/m2),Int_ele_nflx(part/cm2-s-sr),ele_eV(eV from >460eV n/e),Int_ion_eflx,Int_ion_nflx(part/cm2-s-sr),ion_eV (electrons 30.18 keV to 32 eV, ions 30.18 keV to 68 eV' ;printf,u_eflux,' Diff Electron Number Flux - 30 keV to 30 ev' ;printf,u_eflux,' Diff Ion Number Flux - 30 keV to 30 ev' ; 10/08 Emery: only do number of seconds read in particle, not the dimension numsec = n_elements(eflx[*,0]) nendsec = fix(totlen) print,format='(a7,9i8,f7.0)',"numsec ",n_elements(eflx[*,0]),numsec,nendsec,begtime,totlen ;07/10 sometimes totlen (= endtime-curtime) becomes negative so get around it like this if (nendsec lt numsec and nendsec gt 0) then begin endsec = float(begtime[3])*3600. + float(begtime[4])*60. + float(nendsec) if (endsec lt 86339) then numsec=nendsec print,format='(3i8,f8.0)',n_elements(eflx[*,0]),numsec,nendsec,endsec endif for i=0,numsec-1 do begin ;for i=0,n_elements(eflx[*,0])-1 do begin eeflx450(i) = 0 enflx450(i) = 0 ieflx450(i) = 0 inflx450(i) = 0 eeflx(i) = 0 ieflx(i) = 0 enflx(i) = 0 inflx(i) = 0 ; energy channels 0-19 (1-20) are: ; 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ;30.18 20.62 14.04 9.58 6.50 4.42 3.05 2.06 1.41 .984 .992 .679 .462 .317 .213 .145 .100 .068 .046 .032 for j=0,19 do begin if isat ge 16 and j eq 10 then goto,SKIPENFLUX eeflx(i) = eeflx(i) + 1000.*elecal1(j)*ecnt(i,j) enflx(i) = enflx(i) + elecal2(j)*ecnt(i,j) if (j le 12) then begin eeflx450(i) = eeflx450(i) + 1000.*elecal1(j)*ecnt(i,j) enflx450(i) = enflx450(i) + elecal2(j)*ecnt(i,j) ieflx450(i) = ieflx450(i) + 1000.*ioncal1(j)*icnt(i,j) inflx450(i) = inflx450(i) + ioncal2(j)*icnt(i,j) endif if j lt 18 then begin ; Since the two lowest energy ion channels are often contaminated by ambient ; ions drawn in by the spacecraft's potential, ignore these channels in ; calculating the ion energy and number fluxes. These channels are not ignored ; when calculating the differential fluxes so that the contaminated data can be ; seen on the color spectograms. ieflx(i) = ieflx(i) + 1000.*ioncal1(j)*icnt(i,j) inflx(i) = inflx(i) + ioncal2(j)*icnt(i,j) endif SKIPENFLUX: endfor ; 10/08 Emery: ; Calculate mlat and mlt every 1 sec from every 1 min ephemeris data i60 = fix(i/60); ;if (i60*60 eq i) then begin ; printf,u_eflux,format='(a9,f7.0,f7.2)',"UTs_mlat ",uts1m(i60),mlat1m(i60) ;endif rsc1 = (i-i60*60)/60. rsc2 = ((i60+1)*60-i)/60. rmlat = rsc2*mlat1m(i60) + rsc1*mlat1m(i60+1) rmlt = rsc2*mlt1m(i60) + rsc1*mlt1m(i60+1) ; Fix if over 0/24 MLT if (abs(mlt1m(i60)-mlt1m(i60+1)) gt 12) then begin mlt60 = mlt1m(i60) mlt61 = mlt1m(i60+1) if (mlt61 gt mlt60) then rmlt=rsc2*(mlt60+24) + rsc1*mlt61 if (mlt60 gt mlt61) then rmlt=rsc2*mlt60 + rsc1*(mlt61+24) if (rmlt ge 24) then rmlt = rmlt - 24. if (i60*60 eq i-1) then print,format='(a6,3f8.2)',"mlt1m ",mlt60,mlt61,rmlt endif rglat = rsc2*glat1m(i60) + rsc1*glat1m(i60+1) rglat110 = rsc2*glat1101m(i60) + rsc1*glat1101m(i60+1) rglon = rsc2*glon1m(i60) + rsc1*glon1m(i60+1) ; Fix if over 0/360 deg if (abs(glon1m(i60)-glon1m(i60+1)) gt 180) then begin lon60 = glon1m(i60) lon61 = glon1m(i60+1) if (lon61 gt lon60) then rglon=rsc2*(lon60+360) + rsc1*lon61 if (lon60 gt lon61) then rglon=rsc2*lon60 + rsc1*(lon61+360) if (rglon ge 360) then rglon = rglon - 360. if (i60*60 eq i-1) then print,format='(a6,3f8.2)',"lon1m ",lon60,lon61,rglon endif rglon110 = rsc2*glon1101m(i60) + rsc1*glon1101m(i60+1) ; Fix if over 0/360 deg if (abs(glon1101m(i60)-glon1101m(i60+1)) gt 180) then begin lon60 = glon1101m(i60) lon61 = glon1101m(i60+1) if (lon61 gt lon60) then rglon110=rsc2*(lon60+360) + rsc1*lon61 if (lon60 gt lon61) then rglon110=rsc2*lon60 + rsc1*(lon61+360) if (rglon110 ge 360) then rglon110 = rglon110 - 360. if (i60*60 eq i-1) then print,format='(a7,3f8.2)',"lon110 ",lon60,lon61,rglon110 endif rmlon110 = rsc2*mlon1101m(i60) + rsc1*mlon1101m(i60+1) ; Fix if over 0/360 deg if (abs(mlon1101m(i60)-mlon1101m(i60+1)) gt 180) then begin lon60 = mlon1101m(i60) lon61 = mlon1101m(i60+1) if (lon61 gt lon60) then rmlon110=rsc2*(lon60+360) + rsc1*lon61 if (lon60 gt lon61) then rmlon110=rsc2*lon60 + rsc1*(lon61+360) if (rmlon110 ge 360) then rmlon110 = rmlon110 - 360. if (i60*60 eq i-1) then print,format='(a7,3f8.2)',"lon110 ",lon60,lon61,rglon110 endif ralt = rsc2*alt1m(i60) + rsc1*alt1m(i60+1) utsr = uts1m(i60) + float(i-i60*60) ; calculate mean energy from >462eV (mean "E-region" energy - term of Emery) ; See Rich et al., [Annal Geophys, 1987, 87/06 A, 527-534] for use in Robinson et al. ; [1987, JGR, 92, 2565-2570] equs to calculate SigP,H from aurora (NOTE: error in ; eq 14 of Rich et al. - do not divide Eo by 2 for SigH - delete '2' in eq 14.) ; Difference is minor for ions (unless include contamination <68eV), but can ; increase for electrons by almost a factor of 2 since there is lots of ; number flux in electrons <462eV, but very little energy flux. Electrons <600eV ; do not contribute much to the E-region conductivity, but in the cusp or LLBL, ; electrons <600eV contribute a lot to the F-region conductivity, especially in ; cusp reconnection where the downward Poynting flux is large [Li et al., ; 2011, JGR, doi:10.1029/2011JA016566; Knipp et al., 2011, GRL, doi:10.1029/2011GL048302] ele_eV = 0. ion_eV = 0. if (enflx450[i] gt 0.1) then ele_eV = eeflx450[i]/enflx450[i] if (inflx450[i] gt 0.1) then ion_eV = ieflx450[i]/inflx450[i] ; calculate mean energy for electrons from >32eV and for ions >68eV (<68eV contaminated) reg_ele_eV = 0. reg_ion_eV = 0. if (enflx[i] gt 0.1) then reg_ele_eV = eeflx[i]/enflx[i] if (inflx[i] gt 0.1) then reg_ion_eV = ieflx[i]/inflx[i] ;if (abs(rmlat) ge minimum_cglat) then begin ;printf,u_eflux,format='(i5,f7.0,8f8.2,6e12.4)',i,utsr,rglat,rglon,ralt,rglat110,rglon110, $ printf,u_eflux,format='(i5,f7.0,8f8.2,6e12.4,2f7.0)',i,utsr,rglat,rglon,ralt,rglat110,rglon110, $ rmlat,rmlon110,rmlt,eeflx[i]*fac,enflx[i],ele_eV,ieflx[i]*fac,inflx[i],ion_eV,reg_ele_eV,reg_ion_eV ;printf,u_eflux,format='(i5,4e12.4)',i,eeflx[i],ieflx[i],enflx[i],inflx[i] ;printf,u_eflux,format='(5x,10f10.0)',eflx[i,0:19] ;printf,u_eflux,format='(5x,10f10.0)',iflx[i,0:19] ;endif endfor close,u_eflux free_lun,u_eflux ; if the differential number flux (eflx and iflx) have not been computed yet, do so now. ; This only occurs when reading an ascii file which contains only the counts. loc = where(eflx ne 0, count) if count eq 0 then begin for i=0,n_elements(eflx(*,0))-1 do begin IF ISAT LE 15 THEN BEGIN ; SSJ4 for j=0,19 do begin eflx(i,j) = ecnt(i,j)*elecal3(j)/1000.0 iflx(i,j) = icnt(i,j)*ioncal3(j)/1000.0 endfor eflx(i,9) = (eflx(i,9)+eflx(i,10))/2. eflx(i,10) = eflx(i,9) if isat ne 15 then begin iflx(i,9) = (iflx(i,9)+iflx(i,10))/2. iflx(i,10) = iflx(i,9) endif ENDIF ELSE BEGIN ; SSJ5 for j=0,9 do begin eflx(i,j) = ecnt(i,j)*elecal3(j)/1000.0 iflx(i,j) = icnt(i,j)*ioncal3(j)/1000.0 endfor for j=11,19 do begin eflx(i,j) = ecnt(i,j)*elecal3(j)/1000.0 iflx(i,j) = icnt(i,j)*ioncal3(j)/1000.0 endfor ENDELSE endfor endif ; calculate the differential energy flux. Now the user can selected either differential ; energy flux or number flux to be their spectrogram. deeflx = fltarr(n_elements(eflx(*,0)), 20) deiflx = fltarr(n_elements(eflx(*,0)), 20) for i=0,n_elements(eflx(*,0))-1 do begin for j=0,19 do begin deeflx(i,j) = eflx(i,j)*elecal0(j)*1000.0 deiflx(i,j) = iflx(i,j)*ioncal0(j)*1000.0 endfor ; Channel 11 is not meaningful for SSJ5. For consistency, it is not used for J4 either deeflx(i,10) = deeflx(i,9) deiflx(i,10) = deiflx(i,9) endfor ; Next, calculate the average energy flux eave = fltarr(n_elements(eeflx)) iave = fltarr(n_elements(ieflx)) poss = where(enflx gt 0.0, count) if count gt 0 then eave(poss) = eeflx(poss)/enflx(poss) poss = where(inflx gt 0.0, count) if count gt 0 then iave(poss) = ieflx(poss)/inflx(poss) ; If any of our values are less that 1, set them to an absurd value, ; which will be dealt with later. poss = where(eeflx lt 1.0, count) if count gt 0 then eeflx(poss) = 1.0e32 poss = where(ieflx lt 1.0, count) if count gt 0 then ieflx(poss) = 1.0e32 poss = where(eave lt 1.0, count) if count gt 0 then eave(poss) = 1.0e32 poss = where(iave lt 1.0, count) if count gt 0 then iave(poss) = 1.0e32 iic = 0 ; Calculate the integrated number and energy flux from the channels from ; 450 ev to 30 keV integrated over the downward hemisphere assuming ; isotropy oeeflx = fltarr(n_elements(eflx(*,0))) oieflx = fltarr(n_elements(eflx(*,0))) oenflx = fltarr(n_elements(eflx(*,0))) oinflx = fltarr(n_elements(eflx(*,0))) for i=0,n_elements(eflx(*,0))-1 do begin oeeflx(i) = 0 oieflx(i) = 0 oenflx(i) = 0 oinflx(i) = 0 for j=0,12 do begin oeeflx(i) = oeeflx(i) + 1000.*elecal1(j)*ecnt(i,j)*!pi oenflx(i) = oenflx(i) + elecal2(j)*ecnt(i,j)*!pi oieflx(i) = oieflx(i) + 1000.*ioncal1(j)*icnt(i,j)*!pi oinflx(i) = oinflx(i) + ioncal2(j)*icnt(i,j)*!pi endfor endfor ; ; ; Next, calculate the average energy flux from only the 450ev to 30 keV channels oeave = fltarr(n_elements(eeflx)) oiave = fltarr(n_elements(ieflx)) poss = where(oenflx gt 0.0, count) if count gt 0 then oeave(poss) = oeeflx(poss)/oenflx(poss) poss = where(oinflx gt 0.0, count) if count gt 0 then oiave(poss) = oieflx(poss)/oinflx(poss) ; If any of our values are less than 1, set them to an absurd value, ; which will be dealt with later. ; Actually, since we ignore all of the energies below 450 eV, any values ; below 450 would be considered absurd, so set values less than 450 to be ; an even more absurd value, which will be dealt with later. poss = where(oeeflx lt 450.0, count) if count gt 0 then oeeflx(poss) = 1.0e32 poss = where(oieflx lt 450.0, count) if count gt 0 then oieflx(poss) = 1.0e32 poss = where(oeave lt 450.0, count) if count gt 0 then oeave(poss) = 1.0e32 poss = where(oiave lt 450.0, count) if count gt 0 then oiave(poss) = 1.0e32 RETURN END ;----------------------------------------------------------------------------- ; ; dm_ies3.pro - This procedure and its auxilliary procedures open and read ; the SSIES3 Driftmeter from the binary file created by reading the ; SSIES3 EDR file. If the EDR file exists, but the binary file does ; not exist, then the EDR file is read and the binary file is created. ; Since there is limited ephemeris information in the EDR file, the ; ephemeris information is read from the EPH file created at AFRL/Hanscom ; as one-orbit files from AFWA are concatonated into one-day files. ; If the EPH file does not exist, the ephemeris arrays are left blank. ; The IDM values in the EDR file are one-second averages. If 2 or 3 ; second averages are needed for the plot, they are created. ; ;----------------------------------------------------------------------------- PRO dm_ies3, edrfilenm, cor, begtime, inter, dmdata, isatin, $ mlto, mlato, glato, glono, mlatsso, alto, minlat, nsq,$ nodmdata, dm_or_sm_datafound ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ; INPUT: * ; EDRFILENM - the path/name of the SSIES3 EDR file. * ; EDR files processed before Oct 06 have names of the form: * ; ies3_edr_year_day_sssss_eeeee_xxxx.dat * ; where, year=year (e.g. 2004); day=day of year (e.g. 335); * ; sssss=sec of day of start of file; * ; eeeee=sec of day of end of file; * ; xxxx=4 digit spacecraft number (e.g. F16=7554) * ; EDR files process during/after Oct 06 have names of the form:* ; PS.CKGWC_SC.U_DI.A_GP.SIES3-Fxx-R99990-B9999090-APGA_AR.GLOBAL_DD.yearmndy_TP.ssssss-eeeeee_DF.EDR ; where xx=spacecraft number (e.g. 16); year=year * ; mn=number of month; dy=day of month; * ; ssssss=hours/min/sec of start of data; * ; eeeeee=hours/min/sec of end of data * ; If the file name ends with 'Z', then the files is compressed.* ; This program cannot read a compressed file. If the file name* ; ends with 'sav', then the file is an IDL binary file made * ; from the ascii EDR file. This program can read the sav file * ; COR - corotation switch, if=0.0 then remove corotation from * ; drift values, else don't * ; begtime - year,month,day,hour,minute,julian day (all floating point) * ; for start of plot * ; INTER - length of plot in minutes * ; ISATIN - SAT NO., float pt * ; MINLAT - user-given value between 0.0 and 80.0. If it is 0.0 * ; then all latitudes are included in the analysis. * ; NSQ - A SIGNAL, if 0, then plot both hemispheres * ; , if 1, then plot northern hemisphere only * ; , if -1, plot the southern hemisphere only * ; * ; OUTPUTS: * ; dmdata - array of 3 by 12000 dm data points = time (sec of day), vert * ; comp., horz. comp.(12000=6 sample/sec X 60 sec/min X 33 min) * ; However, IES3 EDR file only has 1 sec avg. * ; MLTO - array of 31 magnetic local time at start of each minute,, f.p. * ; MLATO - array of 31 magnetic latitude at start of each minute,, f.p. * ; GLATO - array of 31 geographic latitude at start of each minute,, f.p. * ; GLONO - array of 31 geographic longitude at start of each minute,, f.p.* ; MLATSSO - array of 31 magnetic latitude, sub satellite at start of each * ; minute,, f.p. * ; ALTO - array of 31 altitude at start of each minute, f.p. * ; NODATA - A SIGNAL, if 0, no eof, but no data match in this interval * ; if 1, hit eof and also found no data * ; dm_or_sm_datafound - A SIGNAL, if 0, there is no data in this minlat interval* ; if > 0 logic found an interval of data to plot * ; * ; AUXILLIARY IDL PROCEDURES * ; rdapgafile - reads the IES3 EDR file (ascii) and write values into an * ; IDL binary save file * ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ; Dummy values are returned if not data in interval can be found nodata = 0 dmdata[0,*] = -1 ; start by setting all points to invalid ; Determine if the input file name ends with'.dat' or '.sav'. ; If '.dat', read and convert to an IDL save file. ; if '.sav', skip the conversion and go start to reading the IDL save file. filenm = edrfilenm len_filename = strlen(filenm) ext_filename = strmid(filenm,len_filename-3,3) print,'IES3 EDR file extension: ',ext_filename if ext_filename eq 'sav' then goto,READEDRSAVEFILE if ext_filename ne 'dat' and ext_filename ne 'EDR' then goto,DONEDMIES3 ; Input file must be an ascii EDR file. Read it and convert it to an IDL save file. rdapgafile,filenm filenm2 = strmid(filenm,0,len_filename-3) ; Strip off the 'dat' or 'EDR' from ; the input file name filenm2 = filenm2 +'sav' ; Create the name of the save file filenm = filenm2 READEDRSAVEFILE: restore,filenm ; Contents of Save File: ; imin,min_of_daysv,vhsv,vvsv,satpotsv, yearsv,monthsv,daysv,isat,ephsv, ; plasmadensv,cklsv, ep_densv,engdatasv, ep_anasv, dm_densv, rpa_anasv if begtime[0] ne yearsv[0] or begtime[1] ne monthsv[0] or begtime[2] ne daysv[0] then begin print,'Requested Date does not match EDR Date' print, 'Requested Date =',begtime[0:2] print, 'EDR Date =',yearsv,monthsv,daysv goto,DONEDMIES3 endif print,'Size of input arrays eq Minute in day: ',size(yearsv) print,'Begin Time = ',begtime ; Find minute in saved array that matches (or exceeds) requested beginning minute. index=-1 startmin = 60*begtime[3]+begtime[4] endmin = startmin+inter for i=0,imin-1 do begin if min_of_daysv[i] lt startmin then goto,LOOKNEXTIES3MIN ; Found a time equal to or greater than requested start time if min_of_daysv[i] ge endmin then goto,DONEDMIES3 ; Found a time also less than than time at end of plot ... ; Else quit index = i goto,STARTLOADINGDM LOOKNEXTIES3MIN: endfor STARTLOADINGDM: ; Valid DM data found in time span for plot dm_or_sm_datafound = 0 print,'Data starts at index/min-of-day = ',index,min_of_daysv[index] npts = 0 nodata = 1 dm_or_sm_datafound = 1 NEXTDMIES3: if min_of_daysv[index] ge endmin then goto,DONEDMIES3 for i=0,59 do begin dmdata[0,npts] = min_of_daysv[index]*60+i dmdata[1,npts] = vvsv[i,index] dmdata[2,npts] = vhsv[i,index] npts = npts + 1 endfor index=index+1 goto,NEXTDMIES3 DONEDMIES3: end ;----------------------------------------------------------------------------- ; ; PRO rdapgafile - read the EDR file created by the SSIES3 data processing program, ; known as APGA. Create an IDL save file containing the data ; needed for further processing. ; ; INPUT ; rdapgafile - path/filename of the EDR (ascii) file. Must end with 'dat' or 'EDR' ; ; OUTPUT ; IDL save file with the same path/filename as the input file except the ; 'dat' extension is replaced with 'sav' ;----------------------------------------------------------------------------- ; PRO rdapgafile,filenm,error_flag ; Array Declaration. (Made at start of procedure for readability.) eph = fltarr(6,3,1480) year = intarr(1440) month = intarr(1440) day = intarr(1440) min_of_day= fltarr(1440) satpot = fltarr(15,1440) plasmaden = fltarr(60,1440) vh = fltarr(60,1440) vv = fltarr(60,1440) ckl = fltarr(20,6,1440) ep_den = fltarr(60,1440) ep_ana = fltarr(5,15,1440) rpa_ana = fltarr(8,15,1440) dm_den = fltarr(60,1440) engdata = fltarr(5,1440) ; Create the 'save file' name sl = strlen(filenm) filenm2 = strmid(filenm,0,sl-4) ; Strip off the '.dat' or 'EDR' from the input file name filenm2 = filenm2 +'.sav' ; Create the name of the save file ; Emery 05/07: COMMENT OUT NEXT 2 LINES ADDED IN JAN07 and fail in idl .r syntax error at COUNT=knt ;Result = FILE_SEARCH(filenm2, COUNT=knt ) ;if knt eq 1 then goto,DONEDONE_READEDR error_flag = -1 on_ioerror,NULL get_lun,a_unit openr,a_unit,filenm line=' ' imin = -1 on_ioerror,DONEREADING_EDR print,'Reading SSIES3 EDR Data File' READNEXTMINUTE: IMIN = IMIN + 1 readf,a_unit,line ; First line is a blank line for i=0,115 do begin ; Search for Line 2 readf,a_unit,line ; Line 2 is a title = "RECORD, EDR OF RECORD, DMSP #, DATE, TIME" if strlen(line) gt 8 then begin firstline=strmid(line,0,6) if firstline eq 'RECORD' then goto,FOUNDLINE2 endif endfor print,'Line 2 seems to be incorrect: IMIN=',IMIN print,' ',line GOTO,DONEDONE_READEDR FOUNDLINE2: ; Line 3: header for 1 minute record readf,a_unit,format='(i4,i3,i4,i6,2i2,2x,2i2)',rec_num,edr_num,isat,iyear,imonth,iday,hour,minute year[imin]=iyear month[imin]=imonth day[imin]=iday min_of_day[imin] = 60L*hour+minute ;if imin eq 0 then print,rec_num,edr_num,isat,iyear,imonth,iday,hour,minute,year[imin],month[imin],day[imin] readf,a_unit,line ; Line 4: title = "EPHEMERIS" for i=0,2 do begin ; Lines 5-7 readf,a_unit,a1,a2,a3,a4,a5,a6 ; read ephemeris for minute plus 0, 20 and 40 sec eph[0,i,imin] = a1 ; Geographic Latitude (deg) eph[1,i,imin] = a2 ; Geographic Longitude (deg, E) eph[2,i,imin] = a3 ; Apex Latitude (deg) eph[3,i,imin] = a4 ; Apex Longitude (deg, E) eph[4,i,imin] = a5 ; Apex Local Time (hours) eph[5,i,imin] = a6 ; Satellite Altitude (km) endfor readf,a_unit,line ; Line 8: title = "SATELLITE POTENTIAL, LAST = SOURCE" readf,a_unit,a1,a2,a3,a4,a5,a6,a7,a8 ; Line 9: sat potential at 4 sec intervals satpot[0,imin] = a1 satpot[1,imin] = a2 satpot[2,imin] = a3 satpot[3,imin] = a4 satpot[4,imin] = a5 satpot[5,imin] = a6 satpot[6,imin] = a7 satpot[7,imin] = a8 readf,a_unit,a1,a2,a3,a4,a5,a6,a7,a8 ; Line 10: sat potential plus source flag satpot[8,imin] = a1 satpot[9,imin] = a2 satpot[10,imin] = a3 satpot[11,imin] = a4 satpot[12,imin] = a5 satpot[13,imin] = a6 satpot[14,imin] = a7 satpot_source = fix(a8) readf,a_unit,line ; Line 11: title = "PRIMARY PLASMA DENSITY, THEN SOURCE" for j=0,9 do begin readf,a_unit,a1,a2,a3,a4,a5,a6 ; Lines 12-21: plasma density for 1 minute plasmaden[0+6*j,imin] = a1 plasmaden[1+6*j,imin] = a2 plasmaden[2+6*j,imin] = a3 plasmaden[3+6*j,imin] = a4 plasmaden[4+6*j,imin] = a5 plasmaden[5+6*j,imin] = a6 endfor readf,a_unit,a1 ; Line 22: plasma density source density_source = a1 readf,a_unit,line ; line 23: title = "HORIZONTAL ION DRIFT VELOCS" for j=0,9 do begin ; Line 24-33: one sec avg of horizontal IDM drifts for 1 minute readf,a_unit,a1,a2,a3,a4,a5,a6 vh[0+6*j,imin] = a1 vh[1+6*j,imin] = a2 vh[2+6*j,imin] = a3 vh[3+6*j,imin] = a4 vh[4+6*j,imin] = a5 vh[5+6*j,imin] = a6 endfor readf,a_unit,line ; Line 34: title = "VERTICAL ION DRIFT VELOCS" for j=0,9 do begin ; Lines 35-44: one sec avg of vertical IDM drifts for 1 minute readf,a_unit,a1,a2,a3,a4,a5,a6 vv[0+6*j,imin] = a1 vv[1+6*j,imin] = a2 vv[2+6*j,imin] = a3 vv[3+6*j,imin] = a4 vv[4+6*j,imin] = a5 vv[5+6*j,imin] = a6 endfor readf,a_unit,line ; Line 45: title = "CKL ANALYSES, THEN SOURCE" for j=0,5 do begin ; Line 46-63: CKL analysis at 10 sec intervals for 1 minute readf,a_unit,a1,a2,a3,a4 ckl[0,j,imin] = a1 ckl[1,j,imin] = a2 ckl[2,j,imin] = a3 ckl[3,j,imin] = a4 readf,a_unit,a1,a2,a3,a4,a5,a6,a7,a8 ckl[4,j,imin] = a1 ckl[5,j,imin] = a2 ckl[6,j,imin] = a3 ckl[7,j,imin] = a4 ckl[8,j,imin] = a5 ckl[9,j,imin] = a6 ckl[10,j,imin] = a7 ckl[11,j,imin] = a8 readf,a_unit,a1,a2,a3,a4,a5,a6,a7,a8 ckl[12,j,imin] = a1 ckl[13,j,imin] = a2 ckl[14,j,imin] = a3 ckl[15,j,imin] = a4 ckl[16,j,imin] = a5 ckl[17,j,imin] = a6 ckl[18,j,imin] = a7 ckl[19,j,imin] = a8 endfor readf,a_unit,ckl_source ; Line 64: Source of CKL analysis readf,a_unit,line ; Line 65: title = "EP SWEEP ANALYSES SETS" or "EP AVERAGE DENSITIES" if strmid(line,3,2) eq 'AV' then begin ; Average Densities from Modes C, D and DS for j=0,9 do begin ; Lines 66-75 readf,a_unit,a1,a2,a3,a4,a5, a6 ; read 60 EP avg densities for 1 minute ep_den[0+j*6,imin] = a1 ep_den[1+j*6,imin] = a2 ep_den[2+j*6,imin] = a3 ep_den[3+j*6,imin] = a4 ep_den[4+j*6,imin] = a5 ep_den[5+j*6,imin] = a6 endfor readf,a_unit,line ; Line 76: title = "EP SWEEP ANALYSES SETS" for j=0,2 do begin ; Lines 77-80: up to 3 EP sweep analysis for 1 minute readf,a_unit,a1,a2,a3,a4,a5 ep_ana[0,j,imin] = a1 ; Time ep_ana[1,j,imin] = a2 ; Density ep_ana[2,j,imin] = a3 ; e- temperature ep_ana[3,j,imin] = a4 ; potential ep_ana[4,j,imin] = a5 ; qualifier flag endfor endif else begin ; Sweep Analysis from Modes A, B, BS and E for j=0,14 do begin readf,a_unit,a1,a2,a3,a4,a5 ; Lines 66-80: EP sweep analysis at 4 second intervals ep_ana[0,j,imin] = a1 ; Time ep_ana[1,j,imin] = a2 ; Density ep_ana[2,j,imin] = a3 ; e- temperature ep_ana[3,j,imin] = a4 ; potential ep_ana[4,j,imin] = a5 ; qualifier flag endfor endelse readf,a_unit,line ; Line 81: title = "EP ANALYSES SOURCE" readf,a_unit,ep_source ; Line 82 readf,a_unit,line ; Line 83: title = "RPA SWEEP ANALYSES SETS, THEN SOURCE" for j=0,14 do begin ; Line 84-98: RPA sweep analysis at 4 sec intervals for 1 minute readf,a_unit,a1,a2,a3,a4,a5,a6,a7,a8 rpa_ana[0,j,imin] = a1 ; Sweep Time rpa_ana[1,j,imin] = a2 ; O+ density rpa_ana[2,j,imin] = a3 ; Light Ion Density rpa_ana[3,j,imin] = a4 ; Light Ion Flag rpa_ana[4,j,imin] = a5 ; Ion (O+) Temperature rpa_ana[5,j,imin] = a6 ; Ram ion drift velocity rpa_ana[6,j,imin] = a7 ; Analysis success indicator rpa_ana[7,j,imin] = a8 ; Total Ion density endfor readf,a_unit,rpa_source ; Line 99 readf,a_unit,line ; Line 100: title = "DM ION DENSITY" for j=0,9 do begin ; Lines 101-110: density from IDM data for 1 minute readf,a_unit,a1,a2,a3,a4,a5,a6 dm_den[j*6,imin] = a1 dm_den[1+j*6,imin] = a2 dm_den[2+j*6,imin] = a3 dm_den[3+j*6,imin] = a4 dm_den[4+j*6,imin] = a5 dm_den[5+j*6,imin] = a6 endfor readf,a_unit,line ; Line 111: title = "ENGINEERING DATA" readf,a_unit,a1,a2,a3,a4,a5,a6 ; Line 112: Read Engineering Units for 1 minute engdata[0,imin] = a2 ; ADC Temp engdata[1,imin] = a3 ; SEP Temp engdata[2,imin] = a4 ; DM offset voltage engdata[3,imin] = a5 ; DM mode flag (IES-2: -1=Undefined,0=Normal,1-8=H+,9=FIBA; IES-3:0=Slow, 1=Normal) engdata[4,imin] = a6 ; EP mode (0-6 : A,B,BS,C,D,DS,E) readf,a_unit,line ; Line 113: title = "FILLER" ;print,imin,year[imin],month[imin],day[imin],hour,minute,' ',line readf,a_unit,a1,a2,a3,a4,a5,a6,a7 ; Read Filler values at end of 1 minute record ; Finished reading 1 minute of APGA output goto,READNEXTMINUTE DONEREADING_EDR: close,a_unit free_lun,a_unit ; Save the data iminm1 = imin - 1 ephsv = fltarr(6,3,imin) yearsv = intarr(imin) monthsv = intarr(imin) daysv = intarr(imin) min_of_daysv = fltarr(imin) satpotsv = fltarr(15,imin) plasmadensv = fltarr(60,imin) vhsv = fltarr(60,imin) vvsv = fltarr(60,imin) cklsv = fltarr(20,6,imin) ep_densv = fltarr(60,imin) ep_anasv = fltarr(5,15,imin) rpa_anasv = fltarr(8,15,imin) dm_densv = fltarr(60,imin) engdatasv = fltarr(5,imin) ephsv[*,*,0:iminm1] = eph[*,*,0:iminm1] yearsv[0:iminm1] = year[0:iminm1] monthsv[0:iminm1] = month[0:iminm1] daysv[0:iminm1] = day[0:iminm1] min_of_daysv[0:iminm1] = min_of_day[0:iminm1] satpotsv[*,0:iminm1] = satpot[*,0:iminm1] plasmadensv[*,0:iminm1] = plasmaden[*,0:iminm1] vhsv[*,0:iminm1] = vh[*,0:iminm1] vvsv[*,0:iminm1] = vv[*,0:iminm1] cklsv[*,*,0:iminm1] = ckl[*,*,0:iminm1] ep_densv[*,0:iminm1] = ep_den[*,0:iminm1] ep_anasv[*,*,0:iminm1] = ep_ana[*,*,0:iminm1] rpa_anasv[*,*,0:iminm1] = rpa_ana[*,*,0:iminm1] dm_densv[*,0:iminm1] = dm_den[*,0:iminm1] engdatasv[*,0:iminm1] = engdata[*,0:iminm1] save,imin,min_of_daysv,vhsv,vvsv,satpotsv,$ yearsv,monthsv,daysv,isat,ephsv,$ plasmadensv,cklsv, ep_densv,engdatasv, ep_anasv,$ dm_densv, rpa_anasv, filename=filenm2 DONEDONE_READEDR: print,'Done Reading SSIES EDR file and making IDL save file' end ;----------------------------------------------------------------------------- ; ; dm - This procedure and its auxilliary procedures open and read ; the SSIES2 Driftmeter binary data file, converts the ephemeris ; and data to floating point values, and calculates 1,2 or 3 ; second averages of the ion drift data. ; ;----------------------------------------------------------------------------- PRO dm, filenm$,cor, begtime, inter, dmdata, isatin, $ mlto, mlato, glato, glono, mlatsso, alto, minlat, nsq,$ nodmdata, dm_or_sm_datafound ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ; INPUT: * ; FILENM$ - the name of the driftmeter data file, in ascii * ; COR - corotation switch, if=0.0 then remove corotation from * ; drift values, else don't * ; begtime - year,month,day,hour,minute,julian day (all floating point) * ; for start of plot * ; INTER - length of plot in minutes * ; ISATIN - SAT NO., float pt * ; MINLAT - user-given value between 0.0 and 80.0. If it is 0.0 * ; then all latitudes are included in the analysis. * ; NSQ - A SIGNAL, if 0, then plot both hemispheres * ; , if 1, then plot northern hemisphere only * ; , if -1, plot the southern hemisphere only * ; * ; OUTPUTS: * ; dmdata - array of 3 by 12000 dm data points = time(sec of day),vert * ; comp., horz. comp.(12000=6 sample/sec X 60 sec/min X 33 min) * ; MLTO - array of 31 magnetic local time at start of each minute,, f.p. * ; MLATO - array of 31 magnetic latitude at start of each minute,, f.p. * ; GLATO - array of 31 geographic latitude at start of each minute,, f.p. * ; GLONO - array of 31 geographic longitude at start of each minute,, f.p.* ; MLATSSO - array of 31 magnetic latitude, sub satellite at start of each * ; minute,, f.p. * ; ALTO - array of 31 altitude at start of each minute, f.p. * ; NODATA - A SIGNAL, if 0, no eof, but no data match in this interval * ; if 1, hit eof and also found no data * ; dm_or_sm_datafound - A SIGNAL, if 0, there is no data in this minlat * ; interval; if>0 logic found an interval of data to * ; plot * ; * ; AUXILLIARY IDL PROCEDURES * ; decode_1min_dm * ; corotfix ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * bdata = bytarr(2292) vvert = fltarr(6,60) vhorz = fltarr(6,60) shkp2 = fltarr(60) svap = fltarr(60) gmt0 = fltarr(31) e0 = fltarr(3,31) nodata = 0 dm_or_sm_datafound = 0 isopen = 0 dmdata[0,*] = -1 ; start by setting all points to invalid ;Open the input data file print,filenm$ get_lun,u openr,u,filenm$,ERROR=err if err ne 0 then begin print,"... Error opening dm file <<",filenm$,">>" goto,done_dm endif adata = assoc(u, bytarr(2292)) isopen = 1 ; set up 'null' values in variables for holding info needed to ; skip records to get to the record wanted. lasthrread = -1 lastmnread = -1 hour = -1 minute = -1 lastrecread = 0 istarthr = fix(begtime(3)) istartmin = fix(begtime(4)) istartminofday = 60*istarthr + istartmin iendhr = istarthr iendmin = istartmin + fix(inter) if iendmin gt 59 then begin iendhr = iendhr + 1 iendmin = iendmin - 60 endif iendminofday = 60*iendhr + iendmin num_min_plot = fix(inter) ; Set up couters icycle = 0 idelta_rec = 9999 num_min_saved = 0 nsec_saved = 0 ; begin processing data. This is the top of a loop. READ_DM_AGAIN: ; Read one data record which contains 1 minutes of ephemeris info ; and dm data if not eof(U) then begin on_ioerror,IOERROR_DMREAD bdata = adata(icycle) ; using assoc function to read file. endif else begin print, "DM End of File encountered" if nodata eq 0. and num_min_saved eq 0 then nodata = 1. goto,done_dm endelse isat = 0 year = 0 jday = 0 lasthrread = hour lastmnread = minute hour = 0 minute = 0 geolat = 0. geolong = 0. maglat = 0. mlt = 0. maglon = 0. glatsol = 0. glonsol = 0. glat110 = 0. glon110 = 0. mlat110 = 0. mlon110 = 0. invlat = 0. alt1 = 0 alt2 = 0 bx = 0. by = 0. bz = 0. ex = 0. ey = 0. ez = 0. ssenpot = 0 svbias = 0 svip = 0 srepel = 0 sifree = 0 decode_1min_dm, bdata, vvert, vhorz, shkp2, svap,isat,year,jday,hour,minute, $ geolat, geolong, maglat, mlt, maglon, glatsol, glonsol,$ glat110, glon110, mlat110, mlon110, invlat,alt1,alt2, $ bx,by,bz,ex,ey,ez,ssenpot,svbias,svip,srepel,sifree minofday = 60*hour + minute ;print,format='("Read",3i4,5i6)',jday,hour,minute,minofday,icycle,isat,idelta_rec,lastrecread ; Check for bad minute of data or 'minute of data before requested start' if isat eq 0 or minofday lt istartminofday then goto, SKIPMINUTE if minofday gt istartminofday then begin if idelta_rec gt 1 and idelta_rec lt 9999 then begin ;There must be a data gap. If the program has been skipping records. Re-set ;record position to last read before this point and read forward without skipping idelta_rec = fix(idelta_rec/2) if idelta_rec lt 1 then idelta_rec = 1 icycle = lastrecread+idelta_rec goto,READ_DM_AGAIN endif idelta_rec = 1 endif ; vaid data and time is equal to or greater than requested start time if num_min_saved eq 0 then begin ; skip data if there is no data within 1/2 orbit of start time if minofday gt istartminofday+50 then goto,done_dm ; skip data if the magnetic latitude is too low. if nsq eq 0. then begin if abs(mlat110) lt minlat then goto, SKIPMINUTE endif else begin if nsq*mlta110 lt minlat then goto, SKIPMINUTE endelse begtime(3) = hour begtime(4) = minute istartminofday = 60*hour + minute iendminofday = istartminofday + fix(inter) endif ; Save 1 minute of data ; Save the Ephemeris Info (This is saved even if there are enough minutes ; of drift data already saved. We need N+1 minutes of ephemeris data for ; the corotation calculation.) print,format='(a10,2i4,2f10.1)',"Saving dm ",hour,minute,mlat110,mlt gmt0(num_min_saved) = 3600.*float(hour)+60.*float(minute) mlto(num_min_saved) = mlt mlato(num_min_saved) = mlat110 glato(num_min_saved) = geolat glono(num_min_saved) = geolong mlatsso(num_min_saved)= maglat alto(num_min_saved) = alt1*1.852 e0(0,num_min_saved) = ex e0(1,num_min_saved) = ey e0(2,num_min_saved) = ez ; Check to see if minute read is beyond the requested range if minofday ge iendminofday then goto,done_dm ; Save the Ion Drift data for n=0,59 do begin for nn=0,5 do begin if vhorz(nn,n) gt -3000. and vvert(nn,n) gt -3000. then begin ; Time dmdata[0,nsec_saved] = 3600.*float(hour)+60.*float(minute)+$ float(n)+float(nn)/6. ; Drift dmdata[2,nsec_saved] = vhorz[nn,n] dmdata[1,nsec_saved] = vvert[nn,n] nsec_saved = nsec_saved + 1 endif endfor endfor dm_or_sm_datafound = 1. NODATA = -1. lastrecread = icycle idelta_rec = 1 icycle = icycle + idelta_rec num_min_saved = num_min_saved + 1 goto, READ_DM_AGAIN SKIPMINUTE: ; we got here either becuase of a bad record or record was earlier than ; the start of the range if isat eq 0 then begin ; bad record, read next one idelta_rec = 1 icycle = icycle + idelta_rec if minofday lt iendminofday-1 then goto,READ_DM_AGAIN endif else begin ; record was earlier than start of range lastrecread = icycle if idelta_rec gt 1 then begin i = fix((istartminofday - minofday)/2) if i lt idelta_rec then idelta_rec = i else idelta_rec = idelta_rec - 1 endif else idelta_rec = 1 if idelta_rec lt 1 then idelta_rec = 1 icycle = icycle + idelta_rec goto,READ_DM_AGAIN endelse ; Done reading data for requested time range. GOTO,DONE_DM IOERROR_DMREAD: ; Becuase of logic to skip records, it is possible to try to read a record that is beyond ; the end-of-file. This causes and IOerror condition. Since error may be solved by ; skipping fewer records, the program needs to clear the IOerror condition flag. This is ; done by closing and re-opening the IO unit. close,U openr,u,filenm$ adata =assoc(U, bytarr(2292)) done_dm: if idelta_rec gt 1 and idelta_rec lt 9999 then begin ; If the logic got to this point while trying to skip forward in reading records, then ; cut the size of the number of records skipped and start reading again from the last ; known good record. idelta_rec = fix(idelta_rec/2) if idelta_rec lt 1 then idelta_rec = 1 icycle = lastrecread+idelta_rec goto,READ_DM_AGAIN endif close,U free_lun,u on_ioerror,NULL ; Remove Corotation from the horizontal ion drift (if requested) if num_min_saved lt 2 then stop if cor eq 0.0 and isopen eq 1 then $ COROTFIX, dmdata, nsec_saved, num_min_saved, gmt0, glato, alto, e0 end ;------------------------------------------------------------------------------ ; PRO decode_1min_dm - read 1 minute of driftmeter data from the SSIES Phase 2 ; output file and convert data from packed binary to ; local floating/integer values. ;------------------------------------------------------------------------------ PRO decode_1min_dm, bdata, vvert, vhorz, shkp2, svap,isat,year,jday,hour,minute,$ geolat, geolong, maglat, mlt, maglon, glatsol, glonsol,$ glat110, glon110, mlat110, mlon110, invlat,alt1,alt2,$ bx,by,bz,ex,ey,ez,ssenpot,svbias,svip,srepel,sifree ; ; bdata - array containing on record of SSIES2 data ; vvert - vertical component of ion drift (m/sec), 6 samples/sec X 60 sec ; vhorz - component of ion drift in horizontal plane and perpendicular to ; spacecraft velocity vector, 6 samples/sec X 60 sec ; shkp2 - if not 511, then LLA/LLB value, if eq 511, then H+ mode flag ; svap - aperture potential ; isat - DMSP satellite number (8-20) ; year - calendar year ; jday - julian day of year ; hour - hour of day ; minute - minute of hour ; geolat - geographic latitude of satellite ; geolong - geographic longitude of satellite ; maglat - corr. geomag. lat. of sub-satellite point ; mlt - magnetic local time of satellite ; maglon - corr. geomag. long. of sub-satellite point ; glatsol - geographic latitude of sub-solar point ; glonsol - geographic longitude of sub-solar point ; Begin De-coding the ephemeris block ; Spacecraft ID = F8 to F20 i1 = bdata(1)-48 i2 = bdata(2)-48 if i1 gt 2 then isat = i1 else isat = 10*i1 + i2 ; Sensor ID if isat eq 0 or (bdata(5) ne 68 and bdata(6) ne 77) then begin print,"Not Driftmeter data" print,bdata(0),bdata(1),bdata(2),bdata(3),bdata(4),bdata(5),bdata(6),bdata(7) print,bdata(8),bdata(9),bdata(10),bdata(11),bdata(12),bdata(13),bdata(14),bdata(15) goto, done endif ; year year = bdata(12) year = 1950 + year ; julian day i1 = bdata(13) i2 = bdata(14) jday = 256*i1+i2 ;hour hour = bdata(15) ;minute minute = bdata(16) i1 = bdata(17) i2 = bdata(18) geolat = 256*i1 + i2 geolat = (geolat/10.)-90 i1 = bdata(19) i2 = bdata(20) geolong = 256*i1 + i2 geolong = (geolong/10.) i1 = bdata(21) i2 = bdata(22) maglat = 256*i1 + i2 maglat = (maglat/10.)-90. i1 = bdata(23) i2 = bdata(24) mlt = 256*i1 + i2 mlt = (mlt/10.) i1 = bdata(25) i2 = bdata(26) maglon = 256*i1 + i2 maglon = (maglon/10.) i1 = bdata(27) i2 = bdata(28) glatsol= 256*i1 + i2 glatsol= (glatsol/10.)-90 i1 = bdata(29) i2 = bdata(30) glonsol= 256*i1 + i2 ; for negative values of glatsol, a bias value was used in SSIES Phase I and ; it was not removed in Phase II if glatsol gt 4095 then glatsol = glatsol - 4095 glonsol= (glonsol/10.) i1 = bdata(31) i2 = bdata(32) glat110= 256*i1 + i2 glat110= (glat110/10.)-90 i1 = bdata(33) i2 = bdata(34) glon110= 256*i1 + i2 glon110= (glon110/10.) i1 = bdata(35) i2 = bdata(36) mlat110= 256*i1 + i2 mlat110= (mlat110/10.)-90 i1 = bdata(37) i2 = bdata(38) mlon110= 256*i1 + i2 mlon110= (mlon110/10.) i1 = bdata(39) i2 = bdata(40) invlat = 256*i1 + i2 invlat = invlat /10. i1 = bdata(41) i2 = bdata(42) alt1 = 256*i1 + i2 i1 = bdata(43) i2 = bdata(44) alt2 = 256*i1 + i2 i1 = bdata(45) i2 = bdata(46) bx = long(i1) bx = 256*bx + long(i2) i2 = bdata(47) bx = 256*bx + long(i2) i2 = bdata(48) bx = 256*bx + long(i2) bx = bx/10. - 70000. i1 = bdata(49) i2 = bdata(50) by = long(i1) by = 256*by + long(i2) i2 = bdata(51) by = 256*by + long(i2) i2 = bdata(52) by = 256*by + long(i2) by = by/10. - 70000 i1 = bdata(53) i2 = bdata(54) bz = long(i1) bz = 256*bz + long(i2) i2 = bdata(55) bz = 256*bz + long(i2) i2 = bdata(56) bz = 256*bz + long(i2) bz = bz/10. bz = bz - 70000. i1 = bdata(57) i2 = bdata(58) ex = long(i1) ex = 256*ex + long(i2) i2 = bdata(59) ex = 256*ex + long(i2) ex = ex/100000. ex = ex - 1. i1 = bdata(60) i2 = bdata(61) ey = long(i1) ey = 256*ey + long(i2) i2 = bdata(62) ey = 256*ey + long(i2) ey = ey/100000. ey = ey - 1. i1 = bdata(63) i2 = bdata(64) ez = long(i1) ez = 256*ez + long(i2) i2 = bdata(65) ez = 256*ez + long(i2) ez = ez/100000. ez = ez - 1. ssenpot = bdata(66) svbias = bdata(67)-10 svip = bdata(68)-3 srepel = bdata(69) sifree = bdata(70) ; Begin unpacking one minute of DM data ndm = bdata(71) for i=0,59 do begin for j=0,5 do begin vvert(j,i) = -999990. vhorz(j,i) = -999990. endfor endfor k = 71 for n=0,ndm-1 do begin ; Begin unpacking one second of DM data sec = bdata(k+1) k = k+1 ii = sec for jj=0,5 do begin ; unpack 1 second of v-vertical i1 = bdata(k+1) i2 = bdata(k+2) LL = long(i1) LL = 256*LL + i2 if LL lt 6001 then vvert(jj,ii) = LL*10 - 3000 k = k+2 endfor for jj=0,5 do begin ; unpack 1 second of v-horizontal i1 = bdata(k+1) i2 = bdata(k+2) LL = long(i1) LL = 256*LL + i2 if LL lt 6001 then vhorz(jj,ii) = LL*10 - 3000 k = k+2 endfor i1 = bdata(k+1) i2 = bdata(k+2) shkp2(ii) = long(i1) shkp2(ii) = 256*shkp2(ii) + i2 k = k+2 i1 = bdata(k+1) i2 = bdata(k+2) LL = long(i1) LL = 256*LL + i2 k = k+10 ; move index 2 bytes for data and 8 bytes for zero fill svap(ii) = float(LL)/100. - 19. endfor DONEREADING: DONE: END ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * PRO COROTFIX, FLOW, nsec_saved, NMINSV, gmt, gglat, alt, ecos ; ; COROTFIX - calculate the component of corotation aligned with the ; driftmeter's horizontal component. Subtract this ; corotation from the measured horizontal drift values. ; ; INPUTS: ; flow - array of 1,2 or 3 sec averages of ion drift data. ; (0,*)=time (sec); (1,*) = uncorrected vertical ; ion drift; (2,*) = horizontal ion drift ; nsec_saved - number of drift values in the flow array ; nminsv - number of elements in ephermeris arrays (1 value per minute) ; gmt - time of ephemeris, G.M.T. sec ; gglat - geographic latitude, deg ; alt - altitude, km ; ecos - unit position vector in ECI Coord. ; ; OUTPUTS: ; flow -- flow velocities after removing corotation from ; the horizontal, cross-track component ; ; Comments: ; FORTRAN version written by Marc Hairston, Univ of Texas at Dallas, 1992 ; as part of the program to compute the electrostatic potential across ; the high latitude region from the DMSP IDM data. ; IDL version adapted from FORTRAN version by Fred Rich, Air Force Research ; Lab., Space Veh. Dir, Hanscom AFB, MA, 1998 ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * CORCOR = fltarr(31) PI = 4.*ATAN(1.) RTD = 180./PI ; OMEGA is the angular velocity (radian per second) of the rotating Earth ; used in calculating the corotation of the ionosphere. OMEGA = PI / (12.*3600.) ; Here the program moves the ECI xyz coordinates of the minute to ; ecx, ecy, and ecz and the ECI xyz of the next minute into fcx, fcy, ; and fcz. There should be N+1 sets of ephemeris values, where N=number ; of minutes of IDM data. for IMIN=0,NMINSV do begin i = imin if i eq nminsv then theta = gglat(imin-1)/rtd else theta = gglat(imin)/rtd if i lt nminsv then begin ; position unit vector for beginning of minute ecx = ecos(0,i) ecy = ecos(1,i) ecz = ecos(2,i) ; if there is missing data, then the unit vector will have garbage ; stored in the array. Must replace it will valid, pseudo vector if( abs(ecx) gt 1.0 or abs(ecy) gt 1.0 or abs(ecz) gt 1.0) then begin ecx = 0.01 ecy = 0.01 ecz = 0.9999 endif ; position unit vector for end of minute fcx = ecos(0,i+1) fcy = ecos(1,i+1) fcz = ecos(2,i+1) if( abs(fcx) gt 1.0 or abs(fcy) gt 1.0 or abs(fcz) gt 1.0 ) then begin fcx = 0.0101 fcy = 0.0101 fcz = .9999 endif ; Here the program takes the position of the spacecraft at the present minute ; and at the next minute and determines the velocity vector between them. ; At first the program finds the rotation angle (about the Earth's axis) from ; the y=0 plane to the spacecraft's current position. It then rotates both the ; current and previous positions about the z-axis (which is colinear with the ; Earth's axis) by that angle. Thus the current position is always set at y=0. ; This is a simplifying procedure in that we can now always assume that the ; corotating plasma is flowing in a +y direction, and the only question to be ; determined is the spatial orientation of spacecraft's velocity vector. ; The velocity vector is determined by taking the difference between the compo- ; nents of the current and prior spacecraft positions, then the total length ; of the velocity vector is determined by squaring the components, summing them ; and taking the square root. ;if i lt nminsv then begin ROTANG=ATAN(FCY,FCX) CX=FCX*COS(ROTANG)+FCY*SIN(ROTANG) CY=FCX*SIN(ROTANG)-FCY*COS(ROTANG) BX=ECX*COS(ROTANG)+ECY*SIN(ROTANG) BY=ECX*SIN(ROTANG)-ECY*COS(ROTANG) DELX=CX-BX DELY=CY-BY DELZ=FCZ-ECZ DELTOT=SQRT(DELX*DELX+DELY*DELY+DELZ*DELZ) endif ; The direct angle between the corotating plasma flow vector and the spacecraft ; velocity vector is calculated by taking the arccosine of the delta-y component ; and the hypotenuse (which is the velocity vector itself). By the convention ; used here, the angle is positive if the spacecraft is moving northward and ; negative if the spacecraft is moving southward (the arccosine function always ; returns a positive value). The corotation correction for the horizontal ; IDM flow is a function of the sine of this angle (here done as CORrection). ; The actual plasma corotation velocity at the spacecraft's current location ; (CORotation VELocity) is calculated by multiplying the spacecaft's distance ; from the center of the Earth times the cosine of the latitude times the ; angular velocity of the Earth's rotation. The final corotation correction ; factor (CORotation CORrection) is the multiple of the corotation velocity ; and the correction. The corotation correction factor is then loaded into an ; array for later use. ; DCOSARG=DELY/DELTOT SIGMA=ACOS(DCOSARG) IF DELZ LT 0. then SIGMA=-1.*SIGMA COR=SIN(SIGMA) CORVEL=(6370.0+alt(i))*COS(THETA)*OMEGA*1.E3 CORCOR(I)= -COR*CORVEL endfor ; Here the program fixes the horizontal velocities for corotation factor. ; There is a single loop which passes through all the data. The time of ; each data point is check to be sure it is between two ephemeris times. ; The two correction factors (CORSTRT and CORSTOP) are defined as being the ; correction fraction at the beginning of the current minute and the ; correction factor at the beginning of the subsequent minute. This allows ; a linearly interpolated correction for each data point to be calculated ; and then applied to the flow velocity and loaded back into the ; FLOW array. k = 0 corstrt = corcor(0) corstop = corcor(1) for i=0,nsec_saved-1 do begin if flow(0,i) gt gmt(k+1) then begin k = k+1 corstrt = corcor(k) corstop = corcor(k+1) endif frac = (flow(0,i) - gmt(k))/(gmt(k+1)-gmt(k)) CORFIX = ((CORSTOP-CORSTRT)*FRAC)+CORSTRT FLOW(2,I) = FLOW(2,I) - CORFIX endfor END ;------------------------------------------------------------------------------- ; ; ssm - This procedure and its auxilliary procedures open and read ; the SSM magnetometer binary data file, converts the ephemeris ; and data to floating point values ;------------------------------------------------------------------------------- PRO ssm,smfilen$,begtime,inter,imulf,ssmdata, $ glat0, glon0,mlatss0,alt0,minlat, $ nsq, nossmdata,dm_or_sm_datafound, arrcglalo ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ; INPUT: * ; SMFILEN$ - the name of the magnetometer data file, in ascii * ; BEGTIME(6)- begtime(0) - year * ; - begtime(1) - month * ; - begtime(2) - day all floating point * ; - begtime(3) - hour * ; - begtime(4) - minute * ; - begtime(5) - julian day * ; INTER - length of plot in minutes * ; IMULF - the number of seconds of J4/J5 data to be averaged for output. * ; This number is typically 1, 2 or 3. * ; MINLAT - user-given value between 0.0 and 80.0. If it is 0.0 * ; then all latitudes are included in the analysis. * ; NSQ - A SIGNAL, if 0, then plot both hemispheres * ; , if 1, then plot northern hemisphere only * ; , if -1, plot the southern hemisphere only * ; * ; OUTPUTS: * ; ssmdata - array of 3 by 12000 ssm data points = time, y-component, * ; z component * ; GLAT0 - array of 31 geographic latitude at start of each minute,, f.p. * ; GLON0 - array of 31 geographic longitude at start of each minute,, f.p.* ; MLATSS0 - array of 31 magnetic latitude, sub satellite at start of each * ; minute,, f.p. * ; ALT0 - array of 31 altitude at start of each minute, f.p. * ; NODMDATA - A SIGNAL, = 0, no eof, but no data match in this interval * ; = 1, hit eof and also found no data * ; dm_or_sm_datafound - A SIGNAL, =0, there is no data in this minlat * ; interval; >0 logic found an interval of data to plot* ; * ; AUXILLIARY IDL PROCEDURES * ; cglalo - Convert geodectic lat./long. to corrected geomagnetic coord. * ; * ; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * adata = lonarr(247) t = fltarr(2) idaystart = 0 sec = -1. npt = 0 num_min = 0 ; Pre-set all time elements in data array to -1 to indicate no data. ssmdata[0,*] = -1. ; Open the SSM data file get_lun,U On_IOError,SSMDONE openr,U,smfilen$ ; Set up the "association" between the fixed length records in the input file ; and the internal data array bdata =assoc(U, lonarr(247)) skiptest = 0 icycle = 0 t[0] = 60.*((begtime[3]*60.)+begtime[4]) READMIN: ; Read one minute of data if not eof(U) then begin adata = bdata(icycle) BYTEORDER,ADATA,/NTOHL ;SWAP BYTES 4X4 endif else begin goto, SSMDONE endelse if idaystart eq 0 then idaystart = adata(1) sec = adata[2] sec = sec/1000.+(adata[1]-idaystart)*86400 ; test for starting time if skiptest eq 0 then begin if sec lt t(0) then begin ; Data is prior to the desired start of the plot. Estimate number of ; minutes between data and start time for plot; then skip forward 1/2 ; that number of records. idelta_rec = fix((t(0) - sec)/120.) if idelta_rec lt 1 then idelta_rec = 1 icycle = icycle + idelta_rec goto, READMIN endif else begin if sec gt t(0)+59 then begin skiptest=1 icycle = 0 goto,READMIN endif endelse endif else begin if sec lt t(0) then begin ; Skipping forward did not work. Try again to get the desired ; starting time by reading the data file one record (minute) at ; a time. If 'sec' is found to be greater than the desired ; starting time 't(0)' then the data at the desired starting ; time is missing. Go ahead and plot the next available record. icycle = icycle + 1 goto, READMIN endif endelse ; Proper starting Time found. Convert Ephemeris values to standard float values. glat = adata(3) glat = glat/100. - 90. glng = adata(4) glng = glng/100. glon0[num_min] = glng ; Convert to corrected geomagnetic coordinates CGLALO, arrcglalo, GLAT, GLNG, CGLAT, CGLON mlatss0[num_min] = cglat ; Test for minimum acceptable magnetic latitude. if abs(cglat) lt minlat then begin icycle = icycle + 1 goto, READMIN endif glat0[num_min] = glat ; Test for requested hemisphere. if nsq ne 0 then begin if glat * nsq lt 0.0 then begin if num_min eq 0 then goto, READMIN else goto,SSMDONE endif endif alt = adata(5) alt = alt/10. alt0[num_min] = alt ; Re-set the starting time to the first time that meets all ; the starting criteria. Re-set the ending time to the new ; start time plus the interval for the plot secofday = float(adata(6))/1000. if num_min eq 0 then begin t(0) = secofday t(1) = t(0) + (INTER-1)*60. + 20. endif else if secofday gt t(1) then goto, SSMDONE ; Save data for one minute minpt = 0 print,format='(a5,4i5)','ssm: ',adata(0),adata(1),fix(sec/3600.),fix( (sec/60.) mod 60. ) for i=0,59 do begin ; SAVE A SECOND OF DATA if ( adata(6+minpt) ge 0 and ( npt lt 1859 ) ) then begin ssmdata[0,npt] = adata(6+minpt) ; Time is in milleseconds in file ssmdata[0,npt] = ssmdata[0,npt]/1000. ; convert to seconds ssmdata[1,npt] = adata(8+minpt) ; Delta-BY ssmdata[2,npt] = adata(9+minpt) ; Detta-BZ npt = npt + 1 endif minpt = minpt + 4 endfor num_min = num_min + 1 if ( sec+60 lt t(1) and ( npt lt 12000 ) ) then begin icycle = icycle + 1 goto, READMIN endif SSMDONE: dm_or_sm_datafound = num_min free_lun,U on_ioerror,NULL end ;----------------------------------------------------------------------- pro CGLALO, a, OLDLAT, OLDLON, CGLAT, CGLON ; ; ; THIS SUBROUTINE IS THE HAKURA/GUSTAFSSON COORDINATE TRANSFORMATION ; OF GEOGRAPHIC LATITUDE AND LONGITUDE TO COORECTED GEOMAGNETIC ; LATITUDE AND LONGITUDE. ; ; SYSTEM, KGO REPORT NO. 694 (APRIL 1969), KIRUNA GEOPHYSICAL ; OBSERVATORY, S-981 01 KIRUNA 1, SWEDEN. (BASED ON IGRF65 MODEL) ; ; THIS ROUTINE HAS BEEN CORRECTED AND UPDATED BY MAURIZIO CANDIDI ; AND HERBERT KROEHL OF USDOC/NOAA/DSD, BOULDER, CO 80303. ; THIS ROUTINE HAS BEEN CORRECTED TO ALLOW CALCULATION OF CGM ; COORDINATES RIGHT UP TO THE GEOGRAPHIC POLES. RKR 80/02/20. ; R. K. ROSICH, USDOC/NTIA/ITS, 79/02/23. ; ; MODIFIED TO USE ARGUMENTS BASED ON THE IGRF90 MODEL, BASED UPON WORK OF ; N. E. PAPITASHVILI AND V.O. PAPITASHVILI IN JOURNAL OF ATMOS. AND ; TERRES. PHYSICS, VOL 54, NO. 11/12, PP 1609-1631, 1992. SEE PAPER ; FOR LIMITATIONS OF COORDINATE SYSTEM. COORDINATE SYSTEM IS DEFINED ; FOR THE SURFACE OF THE GEOCENTRIC SPHERE WHICH IS SLIGHTLY DIFFERENT ; THAN ZERO LATITUDE IN GEODEDIC COORDINATES. ; ; F.J.RICH, A.F.PHILLIPS LAB./GEOPHYS.DIR.,BEDFORD,MA 01731. DEC 93 ; ; 12/6/95 ; updated by Radex, Inc., 3 Preston Court, Bedford, MA for F. J. Rich ; using the IGRF95 coefficients and a locally created algorithm to ; reproduce the results of the above paper when the IGRF90 coefficients ; are used. ; ; 4/98 ; Converted from FORTRAN to IDL. Data array made external so routine is ; independent of which version of IGRF is used. ; ;----------------------------------------------------------------------- ; ; ARGUMENTS: ; INPUT: ; A = TABLE OF CORRECTED GEOMAGNETIC LATITUDE AND LONGITUDE PAIRS ; OLDLAT = GEOGRAPHIC LATITUDE (DEGREES, +=N, -=S), ; OLDLON = GEOGRAPHIC LONGITUDE (DEGREES E. LONG.), ; OUTPUT: ; CGLAT = CORRECTED GEOMAGNETIC LATITUDE (DEGREES, +=N, -=S), ; CGLON = CORRECTED GEOMAGNETIC LONGITUDE (DEGREES E. LONG.). ; ;----------------------------------------------------------------------- ; ; corrected magnetic coordinates of geographic north and south poles CGMNLA = a(6408) CGMNLO = a(6409) CGMSLA = a(6410) CGMSLO = a(6411) ;****************************************************************** ; SAVE ORIGINAL INPUTS ;****************************************************************** GGLAT=OLDLAT GGLON=OLDLON ;****************************************************************** ; REDUCE GGLAT TO THE RANGE -90 TO +90 DEGREES ;****************************************************************** point1: if gglat gt 90. then begin gglat = 180. - gglat gglon = gglon+180. goto,point1 endif if gglat lt -90. then begin gglat = -180. - gglat gglon = gglon+180. goto,point1 endif ;*********************************************************************** ; REDUCE GGLON TO .LT. 360 DEGREES (POSITIVE OR NEGATIVE) ;*********************************************************************** GGLON = GGLON mod 360. ;*********************************************************************** ; CHECK TO SEE IF GGLON .GT.0, IF NOT CORRECT TO BE .GT. ZERO ;*********************************************************************** IF GGLON LT 0.0 then GGLON = GGLON + 360.0 I = FIX((90.0 - GGLAT)/2.0) J = FIX(GGLON/10.) + 1 ;*********************************************************************** ; CHECK TO SEE IF I IS WITHIN BOUNDS ;*********************************************************************** IF I GT 0 AND I LT 89 then GOTO,point3 ;*********************************************************************** ; SET TO INTERPOLATE TO THE POLES (N OR S GEOGRAPHIC) ;*********************************************************************** IF J GT 36 then J=J-36 IF J LT 1 then J=J+36 K=J+1 IF K GT 36 then K=K-36 IF I le 0 then begin I=0 X1LAT=CGMNLA X1LON=CGMNLO X2LAT=A( 0+(J-1)*2+(I*72) ) X2LON=A( 1+(J-1)*2+(I*72) ) X3LAT=A( 0+(K-1)*2+(I*72) ) X3LON=A( 1+(K-1)*2+(I*72) ) X4LAT=CGMNLA X4LON=CGMNLO GOTO,point5 endif else begin I=89 X1LAT=A( 0+(J-1)*2+(I-1)*72 ) X1LON=A( 1+(J-1*2)+(I-1)*72 ) X2LAT=CGMSLA X2LON=CGMSLO X3LAT=CGMSLA X3LON=CGMSLO X4LAT=A( 0+(K-1)*2+(I-1)*72 ) X4LON=A( 1+(K-1)*2+(I-1)*72 ) GOTO,point5 endelse ;*********************************************************************** ; STATEMENT 3... CHECK TO SEE IF J IS WITHIN BOUNDS. IF OUT OF ; BOUNDS, RESTORE TO IN BOUNDS ;*********************************************************************** point3: IF J GT 36 then J = J - 36 IF J LT 1 then J = J + 36 ;*********************************************************************** ; FIND CGLAT AND CGLON OF THE FOUR CORNERS OF RECTENGLE WITHIN WHICH ; AND GGLON FALL. THE CORNERS OF THE RECTANGLE ARE X1, X2, X3, AND ; REFERS TO THE NORTHWESTERLY MOST POINT OF THE RECTANGLE. THE LABE ; IN A COUNTERCLOCKWISE MANNER. ;*********************************************************************** X1LAT = A(0+(J-1)*2+(I-1)*72) X1LON = A(1+(J-1)*2+(I-1)*72) X2LAT = A(0+(J-1)*2+I*72) X2LON = A(1+(J-1)*2+I*72) ;*********************************************************************** ; CHECK TO SEE IF THE EASTERLY EDGE OF THE RECTANGLE IS WITHIN THE L ; ARGUMENT BOUNDS OF ARRAY A. ;*********************************************************************** K = J + 1 IF K GT 36 then K = K - 36 X3LAT = A(0+(K-1)*2+I*72) X3LON = A(1+(K-1)*2+I*72) X4LAT = A(0+(K-1)*2+(I-1)*72) X4LON = A(1+(K-1)*2+(I-1)*72) ;*********************************************************************** ; ELIMINATE THE DISCONTINUITIES IN THE X1, X2, X3, AND X4 LONGITUDES ; WOULD OCCUR IF THE PRIME MERIDIAN OF THE CORRECTED GEOMAGNETIC COO ; SYSTEM PASSED THROUGH RECTANGLE X1, X2, X3, X4. ;*********************************************************************** point5: ARR = [X1LON,X2LON,X3LON,X4LON] Y = MAX(ARR) IF ABS(Y-X1LON) GT 180.0 THEN X1LON = X1LON + 360.0 IF ABS(Y-X2LON) GT 180.0 THEN X2LON = X2LON + 360.0 IF ABS(Y-X3LON) GT 180.0 THEN X3LON = X3LON + 360.0 IF ABS(Y-X4LON) GT 180.0 THEN X4LON = X4LON + 360.0 ;*********************************************************************** ; CALCULATE THE SOUTHERN AND THE WESTWARD BOUNDARIES (GEOGRAPHIC LAT ; AND LONGITUDE) OF THE RECTANGLE. ;*********************************************************************** EGGLA = FLOAT(-2*I + 88) WGGLO = FLOAT(10*J - 10) ;*********************************************************************** ; CALCULATE THE COEFFICIENTS OF THE INTERPOLATION EQUATION. THESE ; COEFFICIENTS REPRESENT THE FRACTIONAL PART OF THE BOX THE POINT GG ; GGLON IS FROM THE SOUTHERN AND THE WESTWARD BOUNDARIES OF THE RECT ;*********************************************************************** ALPHA = (GGLAT - EGGLA)/2.0 BETA = (GGLON - WGGLO)/10. ;*********************************************************************** ; NOW CALCULATE THE CORRECTED GEOMAGNETIC COORDINATES THAT CORRESPON ; GEOGRAPHIC COORDINATES OF THE INPUT POINT (GGLAT,GGLON) ;*********************************************************************** CGLAT = X2LAT + ALPHA*(X1LAT - X2LAT) + BETA*(X3LAT - X2LAT) + $ (ALPHA*BETA)*(X4LAT + X2LAT - X1LAT - X3LAT) CGLON = X2LON + ALPHA*(X1LON - X2LON) + BETA*(X3LON - X2LON) + $ (ALPHA*BETA)*(X4LON + X2LON - X1LON - X3LON) ;*********************************************************************** ; CHECK TO DETERMINE IF THE CORRECTED GEOMAGNETIC LONGITUDE IS LESS ; DEGREES. IF IT IS NOT REDUCE THE CALCULATED VALUE BY 360.0 DEGREE ;*********************************************************************** IF CGLON GT 360.0 THEN CGLON = CGLON - 360.0 RETURN END ;----------------------------------------------------------------------- ; PRO SM_NI - procedure to the read the ion density values from the SSIES ; SM file. ; ; INPUT: ; smfile - the name of the SSIES/SM data file, in ascii * ; BEGTIME(6)- begtime(0) - year ; - begtime(1) - month ; - begtime(2) - day all floating point ; - begtime(3) - hour ; - begtime(4) - minute ; - begtime(5) - julian day ; INTER - length of plot in minutes ; ; OUTPUT: ; smnidata - (Time,LogNi) X 1 sec avg X 60 sec/min X 31 min [2,1860] ; ; AUXILIARY PROCEDURES ; decode_1min_sm ;----------------------------------------------------------------------- pro sm_ni,smfile,begtime,inter,smnidata print,'Start sm_ni: ',smfile ; Arrays for one minute of SM data smpower = fltarr(9,60) smdensity = fltarr(60) smrms = fltarr(60) bdata = bytarr(1992) ; Pre-set all time elements in data array to -1 to indicate no data. smnidata[0,*] = -1. ;Open the input data file on_ioerror,END_SMNI get_lun,usmni openr,usmni,smfile adata = assoc(usmni, bytarr(1992)) ; set up 'null' values in variables for holding info needed to ; skip records to get to the record wanted. lasthrread = -1 lastmnread = -1 hour = -1 minute = -1 lastrecread = 0 print,'SM Ni data file opened' ; Enter the starting hour and minute istarthr = begtime[3] istartmin = begtime[4] istartminofday = 60*istarthr + istartmin ; Enter the ending hour and minute iendhr = istarthr iendmin = istartmin + inter if iendmin ge 60 then begin iendhr = iendhr +1 iendmin = iendmin - 60. endif iendminofday = 60*iendhr + iendmin iendminofday = 60*iendhr + iendmin ; Set up counters icycle = 0 idelta_rec = 9999 num_min_saved_for_plot = 0 nsec2plot = 0 ; begin processing data. This is the top of a loop. READAGAIN_SMNI: ; Read one data record which contains 1 minutes of ephemeris info ; and SM data bdata = adata(icycle) ; using assoc function to read file. isat = 0 year = 0 jday = 0 hour = 0 minute = 0 geolat = 0. geolong = 0. maglat = 0. mlt = 0. maglon = 0. glatsol = 0. glonsol = 0. glat110 = 0. glon110 = 0. mlat110 = 0. mlon110 = 0. invlat = 0. alt1 = 0 alt2 = 0 bx = 0. by = 0. bz = 0. ex = 0. ey = 0. ez = 0. ssenpot = 0 svbias = 0 svip = 0 srepel = 0 sifree = 0 decode_1min_sm, bdata, nsm, smpower, smdensity, smrms, $ isat,year,jday,hour,minute, geolat, geolong, maglat, $ mlt, maglon, glatsol, glonsol, glat110, glon110, $ mlat110, mlon110, invlat,alt1,alt2, bx,by,bz,ex,ey,ez,$ ssenpot,svbias,svip,srepel,sifree if isat lt 0 then goto,END_SMNI minofday = 60*hour + minute ; Check for bad minute of data or minute of data before start of plot range if isat eq 0 or minofday lt istartminofday then begin icycle = icycle + 1 ; increment counter without saving minute of data goto, READAGAIN_SMNI endif ; Check for minute of data beyond the plot range if minofday ge iendminofday then goto, END_SMNI ; Minute of data is within plot range. Save it. print,'SM Ni: ',isat,year,hour,minute,invlat im = minofday - istartminofday for ii=0,59 do begin kk = 60*im+ii if smdensity[ii] gt 0. then begin smnidata[0,kk] = 60.*minofday+ii smnidata[1,kk] = smdensity[ii] endif endfor icycle=icycle+1 goto,READAGAIN_SMNI END_SMNI: free_lun,usmni on_ioerror,NULL return end ;=========================================================================== ; ; decode_1min_sm - procedure to convert packed bin into oridinary variables for ; 1 minute of SM data ; ;=========================================================================== PRO decode_1min_sm, bdata, nsm, smpower, smdensity, smrms,isat,$ year,jday,hour,minute, geolat, geolong, maglat,$ mlt, maglon, glatsol, glonsol, glat110, glon110, mlat110, $ mlon110, invlat,alt1,alt2, bx,by,bz,ex,ey,ez,ssenpot,$ svbias,svip,srepel,sifree ; ; nsm - number of valid seconds of data in 1 minute record ; smpower - output from SM comb filters, 9 samples/sec X 60 sec ; smdensity - one second average of 24 samples/sec from SM electrometer, ; Alog10(density), 1 sample/sec X 60 sec ; smrms - variance of one sec avg of alog10(density) ; isat - DMSP satellite number (8-20) ; year - calendar year ; jday - julian day of year ; hour - hour of day ; minute - minute of hour ; geolat - geographic latitude of satellite ; geolong - geographic longitude of satellite ; maglat - corr. geomag. lat. of sub-satellite point ; mlt - magnetic local time of satellite ; maglon - corr. geomag. long. of sub-satellite point ; glatsol - geographic latitude of sub-solar point ; glonsol - geographic longitude of sub-solar point ; Begin De-coding the ephemeris block nsm = 0 ; Spacecraft ID = F8 to F20 i1 = bdata(1)-48 i2 = bdata(2)-48 if i1 gt 2 then isat = i1 else isat = 10*i1 + i2 ; Sensor ID if isat eq 0 or (bdata(5) ne 68 and bdata(6) ne 77) then begin print,"Not SM data" print,bdata(0),bdata(1),bdata(2),bdata(3),bdata(4),bdata(5),bdata(6),bdata(7) print,bdata(8),bdata(9),bdata(10),bdata(11),bdata(12),bdata(13),bdata(14),bdata(15) isat = -1 goto, DONE_NOT_SM_DATA endif ; year year = bdata(12) year = 1950 + year ; julian day i1 = bdata(13) i2 = bdata(14) jday = 256*i1+i2 ;hour hour = bdata(15) ;minute minute = bdata(16) i1 = bdata(17) i2 = bdata(18) geolat = 256*i1 + i2 geolat = (geolat/10.)-90 i1 = bdata(19) i2 = bdata(20) geolong = 256*i1 + i2 geolong = (geolong/10.) i1 = bdata(21) i2 = bdata(22) maglat = 256*i1 + i2 maglat = (maglat/10.)-90. i1 = bdata(23) i2 = bdata(24) mlt = 256*i1 + i2 mlt = (mlt/10.) i1 = bdata(25) i2 = bdata(26) maglon = 256*i1 + i2 maglon = (maglon/10.) i1 = bdata(27) i2 = bdata(28) glatsol= 256*i1 + i2 ; for negative values of glatsol, a bias value was used in SSIES Phase I and ; it was not removed in Phase II if glatsol gt 4095 then glatsol = glatsol - 4095 glatsol= (glatsol/10.)-90 i1 = bdata(29) i2 = bdata(30) glonsol= 256*i1 + i2 glonsol= (glonsol/10.) i1 = bdata(31) i2 = bdata(32) glat110= 256*i1 + i2 glat110= (glat110/10.)-90 i1 = bdata(33) i2 = bdata(34) glon110= 256*i1 + i2 glon110= (glon110/10.) i1 = bdata(35) i2 = bdata(36) mlat110= 256*i1 + i2 mlat110= (mlat110/10.)-90 i1 = bdata(37) i2 = bdata(38) mlon110= 256*i1 + i2 mlon110= (mlon110/10.) i1 = bdata(39) i2 = bdata(40) invlat = 256*i1 + i2 invlat = invlat /10. i1 = bdata(41) i2 = bdata(42) alt1 = 256*i1 + i2 i1 = bdata(43) i2 = bdata(44) alt2 = 256*i1 + i2 i1 = bdata(45) i2 = bdata(46) bx = long(i1) bx = 256*bx + long(i2) i2 = bdata(47) bx = 256*bx + long(i2) i2 = bdata(48) bx = 256*bx + long(i2) bx = bx/10. - 70000. i1 = bdata(49) i2 = bdata(50) by = long(i1) by = 256*by + long(i2) i2 = bdata(51) by = 256*by + long(i2) i2 = bdata(52) by = 256*by + long(i2) by = by/10. - 70000 i1 = bdata(53) i2 = bdata(54) bz = long(i1) bz = 256*bz + long(i2) i2 = bdata(55) bz = 256*bz + long(i2) i2 = bdata(56) bz = 256*bz + long(i2) bz = bz/10. bz = bz - 70000. i1 = bdata(57) i2 = bdata(58) ex = long(i1) ex = 256*ex + long(i2) i2 = bdata(59) ex = 256*ex + long(i2) ex = ex/100000. ex = ex - 1. i1 = bdata(60) i2 = bdata(61) ey = long(i1) ey = 256*ey + long(i2) i2 = bdata(62) ey = 256*ey + long(i2) ey = ey/100000. ey = ey - 1. i1 = bdata(63) i2 = bdata(64) ez = long(i1) ez = 256*ez + long(i2) i2 = bdata(65) ez = 256*ez + long(i2) ez = ez/100000. ez = ez - 1. ssenpot = bdata(66) svbias = bdata(67)-10 svip = bdata(68)-3 srepel = bdata(69) sifree = bdata(70) ; Begin unpacking one minute of SM data nsm = bdata(71) smdensity[*] = -999990. k = 71 for n=0,nsm-1 do begin ; Begin unpacking one second of SM data sec = bdata(k+1) k = k+1 ii = sec ;for jj=0,8 do begin ; unpack 1 second of output of sm comb filters ; i1 = bdata(k+1) ; i2 = bdata(k+2) ; LL = long(i1) ; LL = 256*LL + i2 ; smpower(jj,ii) = float(LL)/1000. ; k = k+2 ;endfor ; since the filter output is usually bad, we will skip processing the ; output of the filters. k = k + 18 ; unpack 1 second average of Log10(ion density [cm^-3]) i1 = bdata(k+1) i2 = bdata(k+2) LL = long(i1) LL = 256*LL + i2 if LL lt 60001 then smdensity[ii] = float(LL)/10000. k = k+2 ; unpack 1 second of rms of 1 sec avg of Log10(density) i1 = bdata(k+1) i2 = bdata(k+2) smrms(ii) = long(i1) smrms(ii) = (256*smrms(ii) + i2)/10000. k = k+2 ; zero file k = k + 9 endfor DONE_NOT_SM_DATA: RETURN END ;----------------------------------------------------------------------- ; PRO SM_NI_IES3 - procedure to the read the ion density values from the SSIES3 ; EDR file. ; ; INPUT: ; edrfilenm - the name of the SSIES3 EDR data file (ascii) or ; the name of the SSIES3 EDR IDL Save file (binary) ; BEGTIME(6)- begtime(0) - year ; - begtime(1) - month ; - begtime(2) - day all floating point ; - begtime(3) - hour ; - begtime(4) - minute ; - begtime(5) - julian day ; INTER - length of plot in minutes ; ; OUTPUT: ; smnidata - (Time,LogNi) X 1 sec avg X 60 sec/min X 31 min [2,1860] ; ; AUXILIARY PROCEDURES ; decode_1min_sm ;----------------------------------------------------------------------- pro sm_ni_ies3,edrfilenm,begtime,inter,smnidata print,'Start sm_ni_ies3: ',edrfilenm ; Pre-set all time elements in data array to -1 to indicate no data. smnidata[0,*] = -1. ; Determine if the input file name ends with'.dat' or '.sav'. ; If '.dat', read and convert to an IDL save file. ; if '.sav', skip the conversion and go start to reading the IDL save file. filenm = edrfilenm len_filename = strlen(filenm) ext_filename = strmid(filenm,len_filename-3,3) ;print,'IES3 EDR file extension: ',ext_filename if ext_filename eq 'sav' then goto,READEDRSAVEFILESM if ext_filename ne 'dat' and ext_filename ne 'EDR' then goto,DONESMIES3 ; Input file must be an ascii EDR file. Read it and convert it to an IDL save file. rdapgafile,filenm filenm2 = strmid(filenm,0,len_filename-3) ; Strip off the 'dat' from the input file name filenm2 = filenm2 +'sav' ; Create the name of the save file filenm = filenm2 READEDRSAVEFILESM: restore,filenm ; Contents of Save File: ; imin,min_of_daysv,vhsv,vvsv,satpotsv, yearsv,monthsv,daysv,isat,ephsv, ; plasmadensv,cklsv, ep_densv,engdatasv, ep_anasv, dm_densv, rpa_anasv if begtime[0] ne yearsv[0] or begtime[1] ne monthsv[0] or begtime[2] ne daysv[0] then begin print,'Requested Date does not match EDR Date' print, 'Requested Date =',begtime[0:2] print, 'EDR Date =',yearsv,monthsv,daysv goto,DONESMIES3 endif print,size(yearsv) print,begtime ; Find minute in saved array that matches (or exceeds) requested beginning minute. index=-1 startmin = 60*begtime[3]+begtime[4] endmin = startmin+inter for i=0,imin-1 do begin if min_of_daysv[i] lt startmin then goto,LOOKNEXTIES3MIN ; Found a time equal to or greater than requested start time if min_of_daysv[i] ge endmin then goto,DONESMIES3 ; Found a time also less than than time at end of plot ... ; Else quit index = i goto,STARTLOADINGSM LOOKNEXTIES3MIN: endfor STARTLOADINGSM: ; Valid DM data found in time span for plot print,index,min_of_daysv[index] npts = 0 NEXTDMIES3: if min_of_daysv[index] ge endmin then goto,DONESMIES3 for i=0,59 do begin smnidata[0,npts] = min_of_daysv[index]*60+i smnidata[1,npts] = alog10(plasmadensv[i,index]) npts = npts + 1 endfor index=index+1 goto,NEXTDMIES3 DONESMIES3: return end ;------------------------------------------------------------------------------ ; ; MAIN PROGRAM ; ;------------------------------------------------------------------------------ close,/all ; Just in case the last IDL excution aborted. ; Initialize variables ; comment out by Lei in 01/07 (Emery in 05/07): ;eflx = fltarr(1800,20) ; 30 minute (max plot length) x 60sec/min x 1 sample/sec, 20 channels ;iflx = fltarr(1800,20) ;ecnt = fltarr(1800,20) ;icnt = fltarr(1800,20) ;dmdata = fltarr(3,12000) ; 12000=(Time,V_Vert,V_horz) X 6 samples/sec x 60 sec/min X 33.3333 min ;ssmdata = fltarr(3,1860) ; (Time,Delta-By,Delta-Bz) X 1 sec avg X 60 sec/min X 31 min ;smnidata = fltarr(2,1860) ; (Time,LogNi) X 1 sec avg X 60 sec/min X 31 min print,' ' print,' Welcome to Cenergy.pro' print,' An IDL program to make publication quality plots of DMSP SSJ/SSIES/SSM data' print,' Version August 2004' ; Ephemeris data variables ; comment out by Lei in 01/07 (Emery in 05/07): ;mlt = fltarr(31) ;mlat = fltarr(31) ;mlt2 = fltarr(31) ;mlat2 = fltarr(31) ;glat = fltarr(31) ;glon = fltarr(31) ;mlatss = fltarr(31) ;alt = fltarr(31) ; random date stuff ; The 28 for dayofmon is changed to 29 if the year is leap. ; monname = 'JanFebMarAprMayJunJulAugSepOctNovDec' dayofmon = [31,28,31,30,31,30,31,31,30,31,30,31] ; blank = [' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '] ; set up Energy Y-Axis yb = 0.0 ye = 419.0 ytick = 2 ytickval = [91.0,223.0,355.0] yticklab = ['10!A2','10!A3','10!A4'] ; set up Electron and Ion color chart yaxis cyb = 0.0 cye = 15.0 ctick = 5 ctickval = fltarr(ctick+1) for i=0,5 do ctickval(i) = cye - (cye*0.95/5.0)*float(5-i) cticke = ['10!E2',' ','10!E4',' ','10!E6',' '] cticki = ['10!E0',' ','10!E2',' ','10!E4',' '] print, ' ' output_type = 0 inputerror=0 ps_encapsulated_option=0 get_file_names, partfile, removebackground, driffile, ssmffile, smfile, psfile,$ datafile_or_keyboard_input, output_type, dm_or_ssm_option,inputerror,$ ps_encapsulated_option if inputerror lt 0 then goto,final1 ; if SSM data is going to be plotted, then we need to load the table of ; values for converting from Geodectic coord. to Corrected Geomagnetic coord. if dm_or_ssm_option eq 2 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin arrcglalo = fltarr(6416) ; 2 X 36 X 89 + 4 for poles + 4 zero-fill ; if output_type gt 0 then filenam = "./cglalo2k.dat" else $ ; filenam =".\cglalo2k.dat" if output_type gt 0 then filenam = "./cglalo2005.dat" else $ filenam =".\cglalo2005.dat" on_ioerror,CGLALONOTFOUND is_cglalo = 0 get_lun,u openr,u,filenam is_cglalo = 1 ; open was successful CGLALONOTFOUND: on_ioerror,NULL if is_cglalo eq 0 then begin ; print,"cglalo95.dat was not found in the default directory" print,"cglalo2k.dat was not found in the default directory" print,"Enter the path/name of the file with values for Geodectic to Corr. Geomag" print," conversion" read,filenam print,filenam openr,u,filenam is_cglalo = 1 ; open was successful endif print,"Please Wait ... Obtaining the CGM Transformation array" for i=0,801 do begin readf,u,a0,a1,a2,a3,a4,a5,a6,a7 arrcglalo[8*i]=a0 arrcglalo[8*i+1]=a1 arrcglalo[8*i+2]=a2 arrcglalo[8*i+3]=a3 arrcglalo[8*i+4]=a4 arrcglalo[8*i+5]=a5 arrcglalo[8*i+6]=a6 arrcglalo[8*i+7]=a7 endfor close,u free_lun,u print,"Finished obtaining CGM array" endif ; If you request a postscript (encapsulated postscript) file, then you ; should have named the output file with a .ps (.eps) file extension. ; Break the filename down into the path/file_name and extension, ; so the program can add a counter onto the each name and the different ; passes will distinguished from each other. if abs(output_type) eq 1 then begin j=-1 for i=0,strlen(psfile)-1 do if strmid(psfile,i,1) eq '.' then j=i if j ge 0 then begin filemain = strmid(psfile,0,j) filext = strmid(psfile,j,strlen(psfile)-j) endif else begin filemain = psfile if ps_encapsulated_option eq 0 then filext = '.ps' else filext = '.eps' endelse endif ; Driftmeter or no driftmeter if (strmid(driffile,0,4) ne 'null') and $ (strmid(driffile,0,4) ne 'NULL') then nodm = 0 else nodm = 1 ; ; Magnetometer or no magnetometer ; if (strmid(ssmffile,0,4) ne 'null') and $ (strmid(ssmffile,0,4) ne 'NULL') then nossm = 0 else nossm = 1 ; SSIES/SM Ni data or no data if (strmid(smfile,0,4) ne 'null') and $ (strmid(smfile,0,4) ne 'NULL') then nosmni = 0 else nosmni = 1 z=0.0 begtime = intarr(6) get_time_int, dayofmon, begtime, endtime, length, satel,cor, ycol, yfil, $ maxvel, max_nt, minlat, nsq, nodm, nossm, vopt,$ datafile_or_keyboard_input, portrait_or_land_option, diswitch, flipi,$ vcount, vlabc, vstr, vlabx, vlines, intsw, imul,output_type if portrait_or_land_option eq 0 then begin ; For landscape the xsize and ysize of the page in inches is as follows: change=1.0 ; if this value is decrease, the size of the plot will be decreased xs = 10.0 ys = 7.5 endif else begin ; For portrait the sizes are as follows: change=1.0 xs = 7.5 ys = 7.5*0.75 ; Since we are keeping the same proportions as landscape (in x and y) ; we have to make sure that the letters are scaled down to the right ; size also, so the scale has to be modified by 0.75 to correct this. ; But we need the original change for the device call, so we store it ; in changep. changep = change change = 0.75*change endelse ; The absolute longest time (ABSLEN) for each plot is LENGTH abslen = length isat = float(satel) imulf = float(imul) ; Compute the current time (the first time to look at: CURTIME) in ; seconds curtime = float(begtime[3])*3600.0 + float(begtime[4])*60.0 ; Keep track of the number of plots to put out. One page for each plot. pagenum = 0 ; If we open and close the J4/J5 file every time we want to search for a ; specific time, the program would take forever to run, since J4/J5 files ; can be as long as people want them to be. So, we open the file once ; and save our place in it, then at the end, we close it. ; if j4op = 0, then the file will be opened and j4op will be set to 1 ; if j4cl = 1, then the file will be closed and no data retrieved j4op = 0.0 j4cl = 0.0 ; This loop is set up so the program will plot out from the start time ; to the end time. But, if the end time falls in the middle of a plot, ; the program will plot the entire plot. So the end time is really for ; the start of a plot (if the end time of a plot is after the end time ; of the program, the program will end.) ; This is done by the programmer's choice. It can be changed relatively ; easily, if one wanted the end time to be the absolute end time. ; ; The program keeps track of the "time" of the file with the variables ; CURTIME and BEGTIME. CURTIME is the time in seconds from the start of ; the day. BEGTIME is an array which contains the date and time (see ; above description.) ; ; 10/08 Emery: change lt to le ;while (curtime lt endtime) do begin while (curtime le endtime) do begin ; add here by JH Lei in 01/07 (Emery in 05/07) just after: while (curtime lt endtime) do begin ; Initialize variables eflx = fltarr(1800,20) ; 30 minute (max plot length) x 60sec/min x 1 sample/sec, 20 channels iflx = fltarr(1800,20) ecnt = fltarr(1800,20) icnt = fltarr(1800,20) dmdata = fltarr(3,12000) ; 12000=(Time,V_Vert,V_horz) X 6 samples/sec x 60 sec/min X 33.3333 min ssmdata = fltarr(3,1860) ; (Time,Delta-By,Delta-Bz) X 1 sec avg X 60 sec/min X 31 min smnidata = fltarr(2,1860) ; (Time,LogNi) X 1 sec avg X 60 sec/min X 31 min ; Ephemeris data variables mlt = fltarr(31) mlat = fltarr(31) mlt2 = fltarr(31) mlat2 = fltarr(31) glat = fltarr(31) glon = fltarr(31) mlatss = fltarr(31) alt = fltarr(31) ; 10/08 Emery: additional ephemeris data variables dglat = fltarr(31) dglon = fltarr(31) dglat110 = fltarr(31) dglon110 = fltarr(31) dmlon110 = fltarr(31) salt = fltarr(31) uts = fltarr(31) ; INTER is set at first( in MAIN) to the absolute length of a plot. ; It will remain this value, unless the user has specified a min lat ; to plot OR the beginning and end time cannot be found in the data ; file which is tested in the subroutine dm_plot. Then the ; procedure which reads the driftmeter data will change the value to the ; true time interval of the plot. If INTER is returned as zero, no ; plot is made. inter = abslen realtime = float(begtime) ; Total time in seconds remaining in the requested time: totlen = endtime - curtime ; ; Logic to set up switches, namely, nodmdata if begtime[4] gt 9 then $ print, 'checking data for '+strcompress(string(begtime(3))+':'+ $ string(begtime(4)),/remove_all) $ else $ print, 'checking data for '+strcompress(string(begtime(3))+':0'+ $ string(begtime(4)),/remove_all) ; Call the procedures to acquire the data to be plotted ; if dm_or_ssm_option eq 3.0 then begin nodmdata = 1.0 inter = 0.0 goto, cont3 endif ; ict = 0.0 nodmdata = 0.0 dm_or_sm_datafound = 0.0 while (dm_or_sm_datafound eq 0.0) do begin ; DM data to be included in the plot (either alone or with SSM & SM data) if dm_or_ssm_option eq 1 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin nodmdata = 0.0 IF satel lt 16 then begin ; Read the SSIES/SSIES2 IDM data file dm, driffile,cor, realtime, inter, dmdata, isat, $ mlt, mlat, glat, glon, mlatss, alt, minlat, nsq,$ nodmdata, dm_or_sm_datafound endif else begin ; Read the SSIES3 IDM data from the EDR file; Read the ephemeris information ; from the EPH file. dm_ies3, driffile,cor, realtime, inter, dmdata, isat, $ mlt, mlat, glat, glon, mlatss, alt, minlat, nsq,$ nodmdata, dm_or_sm_datafound endelse endif ; SSM data to be included in the plot (either alone or with DM & SM data) if dm_or_ssm_option eq 2 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin nodmdata = 0.0 ; Read the SSM data file ssm,ssmffile,realtime,inter,imulf,ssmdata, $ glat, glon,mlatss,alt, $ minlat,nsq, nodmdata,dm_or_sm_datafound, $ arrcglalo endif ; SM Ni data to be included in the plot if dm_or_ssm_option eq 4.0 then begin IF satel lt 16 then begin ; Read the SSIES/SM data file sm_ni,smfile,realtime,inter,smnidata endif else begin ; Read the SSIES3 SM data from the EDR file sm_ni_ies3,smfile,realtime,inter,smnidata endelse endif ; begtime = fix(realtime) ict = ict + 1.0 endwhile ; ; If no data is returned, no plot is made ; We are able to test dmdata because either driftmeter or magnetometer ; is stored there during job execution, never both <<<<<<< EXCEPT FOR OPTION 4! ; cont3: poss = where(dmdata(0,*) ne 0.0, count) if count eq 0 then inter = 0.0 ; If there was a time interval, the beginning time was returned as a ; float value in REALTIME. We convert it back to integer. if inter ne 0.0 then begtime = fix(realtime) ; Want to make sure that the longest time period plotted is ABSLEN, ; and set the amount of time we want to plot to LENGTH if (abslen gt inter) and (inter ne 0.0) then $ length = inter $ else length = abslen if (inter eq 0.0) and (nodmdata eq 1.0) then $ print, 'No Driftmeter data found in this interval.' $ else $ if (inter eq 0.0) and (nodmdata eq 0.0) then $ print, 'No data in the Latitude Range Selected.' $ else begin if nodm eq 0 then print,'Data for Driftmeter found.' if nossm eq 0 then print,'Data for Magnetometer found.' endelse ; Call the procedure for the J4/J5 data if we have an actual ; time interval dumlat = minlat dumnsq = nsq dumlen = totlen/60. inter = length plotparsw = 0 ; ; Attempt to iterate thru inter no. of minutes of data until ; the new switch plotparsw is set higher than 0 from icount in ; subroutine particle. It should take a group of "inter" minutes ; at a time. ; if (strmid(partfile,0,4) eq 'null') or $ (strmid(partfile,0,4) eq 'NULL') $ then goto, cont4 ; ; if inter ne 0.0 then begin ecnt = fltarr(1800,20) ecnt[0:1799,0:19] = 0 icnt[0:1799,0:19] = 0 iyear = fix(begtime[0]) ; Read data from the SSJ4/SSJ5 file particle, partfile, realtime, inter, imulf, eflx, ecnt, iflx, $ icnt, isat, mlat2, mlt2, j4op, j4cl, $ ; 10/08 Emery: add more ephemeris data ; dumlat, dumnsq, dumlen, plotparsw, removebackground, iyear dumlat, dumnsq, dumlen, plotparsw, removebackground, iyear, $ dglat,dglon,salt,dglat110,dglon110,dmlon110,uts ; 10/08 Emery: add specific ephemeris data from 0000 UT of the next day file for calc_flx iendminofd = fix(begtime[3])*60. + fix(begtime[4]) + fix(abslen) if (iendminofd ge 1439) then begin ; Set values of 31st uts etc when endtime 2359 (1439 min) - want values at 2400 UT on next file uts[30] = 86400. ; This is a 'hack' to get the ephemeris data correct for the last point in the file. ; It uses a correction file (ex: cenergy_hack_F15_2005_05_15.txt) containing the first 2 minutes of SSJ data for the following day, ; for a particular satellite. The naming convention is the 'hack' file is named for the day it is used to correct, NOT the day the data comes from. s = string(format='(a14,i2.2,a1,i4.4,i2.2,i2.2,a4)','cenergy_hack_f',isat,'_',begtime[0],begtime[1],begtime[2],'.txt') ; June 2012 Emery: Cannot find where these prints go (other prints are in docenergy2.sh) if FILE_TEST(s) then begin GET_LUN, U OPENR, U, s header = STRARR(2) READF, U, header number = 1 uttime = 1 ; data = dblarr(16) data = fltarr(16) st = '' ; READF, U,number, uttime,st, data,FORMAT='(i0,i0,A1,8f0,:)' READF, U,number, uttime,st, data,FORMAT='(i0,i0,A1,f0,:)' print, data, number, uttime, st dglat[30] = data(0) dglon[30] =data(1) salt[30]=data(2) dglat110[30]=data(3) dglon110[30]=data(4) mlat2[30]=data(5) dmlon110[30]=data(6) mlt2[30]=data(7) print,'Correction applied' print,s CLOSE, U FREE_LUN, U endif else begin print, 'Correction file not found, duplicating previous point for last data point' print,s dglat[30]=dglat[29] dglon[30]=dglon[29] salt[30]=salt[29] dglat110[30]=dglat110[29] dglon110[30]=dglon110[29] mlat2[30]=mlat2[29] dmlon110[30]=dmlon110[29] mlt2[30]=mlt2[29] endelse ; values for j4f1505136: ; dglat[30] = 38.5 ; dglon[30] = 139.8 ; salt[30] = 454. ; dglat110[30] = 42.7 ; dglon110[30] = 139.1 ; dmlon110[30] = 210.6 ; mlat2[30] = 35.8 ; mlt2[30] = 9.626 endif ; 10/08 Emery: TEMP print print,FORMAT='(a6,31f7.0)',"utsec ",uts print,FORMAT='(a5,31f7.2)',"mlat ",mlat2 endif begtime = fix(realtime) ; cont4: ; if inter ne 0.0 then begin ; Compute where the tick marks should be and what their values should be begtime = fix(realtime) length = inter if dm_or_ssm_option eq 1.0 then begin if (nodm eq 1.0) or (nodmdata eq 1.0) then begin ; Switches indicate no dm data in this interval. Zero out data locations. for isec=0, inter*60-1 do $ dmdata(0,isec)=float(isec+begtime(3)*3600.0+begtime(4)*60.0) dmdata(0,inter*60:n_elements(dmdata(0,*))-1) = 0.0 endif endif else if dm_or_ssm_option eq 2.0 then begin if (nossm eq 1.0) or (nodmdata eq 1.0) then begin ; Switches indicate no ssm data in this interval. Zero out data locations. ; for isec=0, inter*60-1 do $ dmdata(0,isec)=float(isec+begtime(3)*3600.0+begtime(4)*60.0) dmdata(0,inter*60:n_elements(dmdata(0,*))-1) = 0.0 endif endif tick_stuff, begtime, length, nticks, tickval, ticklab, bt, et, imul ; The J4/J5 data exists then go ahead and plot it all out if (count ne 0) or (plotparsw ne 0) then begin ; Want to limit the data to the amount of time selected, and not include ; any values which may have lingered around, so cut the arrays down to ; the number of relevant points npts = fix(inter) tot = npts*60-1 eflx2 = eflx(0:tot,*) iflx2 = iflx(0:tot,*) ecnt2 = ecnt(0:tot,*) icnt2 = icnt(0:tot,*) if plotparsw ne 0 then begin ; Interpolate the data and calculate the flux and energy ; print,"Beginning conversion of J4/J5 data from counts to flux" JGF_LABEL = ' ' calc_flx, ecnt2, icnt2, isat, eeflx, ieflx, eave, iave, eflx2, iflx2, $ ; Emery 05/07: added begtime to calc_flx ; Emery 10/08: added totlen, dumlat and ephemeris data to calc_flx oeeflx, oieflx, oenflx, oinflx, deeflx, deiflx, output_type, iyear, $ ; begtime, JGF_LABEL begtime,JGF_LABEL,totlen,dumlat,uts,dglat,dglon,salt,dglat110, $ dglon110,dmlon110,mlat2,mlt2 ; print,'main: calc_flx: JGF_LABEL = ', JGF_LABEL ; print,"End counts to flux conversion" ; If user selects to display differential energy flux instead of differential ; Number flux, change arrays, so the colors come out right. if diswitch eq 1 then begin eflx2 = deeflx iflx2 = deiflx endif print,"Start creation of bitmap for spectograms" interpo, eflx2, iflx2, npts, ele2d, ion2d, diswitch,output_type, $ dm_or_ssm_option print,"End creation of bitmap for spectograms" endif else begin print,"Particle Data Not Found" endelse ; Increment the page number pagenum = pagenum + 1 ; If there is a postscript file, and more than one page, edit the file ; name to include the page counter if (psfile ne '') and (pagenum gt 1) then $ psfile=filemain+strcompress(string(pagenum),/remove_all)+filext ; These are used for the adding of lines and points, if no ps file is ; entered postps = '' first_time = 1 if psfile ne '' then postps = psfile ; !P.THICK = 1.0 ; Standard Line Weight if abs(output_type) eq 1 or abs(output_type) eq 2 then begin ; The intent is to make a post_script file or HP Laserjet print file ; ; This sets the font to be hardware drawn for post-script !P.Font = 0 ; Graphics device specific font to be used ; 07/10 Liam Kilcommons to disable plotting: GOTO, loop_it ; Disable plotting; just skip to the next step in the while loop if abs(output_type) eq 1 then begin ; post-script or encapsulated post-script set_plot, 'ps', /copy, /interpolate if ps_encapsulated_option eq 0 then begin ; standard post-script file ; color or no color if ycol eq 0.0 then begin ; no color if portrait_or_land_option eq 0 then begin ; landscape device, file=postps, /landscape, bits=8, xsize=xs*change, $ ysize=ys*change, yoff=11.0-(11.0-xs*change)/2.0, $ xoff=(8.5-ys*change)/2.0, /inches, /times, encapsulated=0, preview=0 endif else begin ; portrait device, file=postps, /portrait, bits=8, xsize=xs*changep, $ ysize=ys*changep, yoff=(11.0-ys*changep)/2.0, $ xoff=(8.5-xs*changep)/2.0,/inches, /times, encapsulated=0, preview=0 endelse endif else begin ; color if portrait_or_land_option eq 0 then begin ; landscape device, file=postps, /landscape, /color, bits=8, $ xsize=xs*change, ysize=ys*change, $ yoff=11.0-(11.0-xs*change)/2.0, xoff=(8.5-ys*change)/2.0, $ /inches, /times, encapsulated=0, preview=0 endif else begin ; portrait device, file=postps, /portrait, bits=8, xsize=xs*changep, $ ysize=ys*changep, yoff=(11.0-ys*changep)/2.0, /color, $ xoff=(8.5-xs*changep)/2.0, /inches, /times, encapsulated=0, preview=0 endelse endelse endif else begin ; encapsulated post-script ; color or no color if ycol eq 0.0 then begin ; no color if portrait_or_land_option eq 0 then begin ; landscape device, file=postps, /landscape, bits=8, xsize=xs*change, $ ysize=ys*change, yoff=11.0-(11.0-xs*change)/2.0, $ xoff=(8.5-ys*change)/2.0, /inches, /times, preview=1, encapsulated=1 endif else begin ; portrait device, file=postps, /portrait, bits=8, xsize=xs*changep, $ ysize=ys*changep, yoff=(11.0-ys*changep)/2.0, $ xoff=(8.5-xs*changep)/2.0,/inches, /times, preview=1, encapsulated=1 endelse endif else begin ; color if portrait_or_land_option eq 0 then begin ; landscape device, file=postps, /landscape, /color, bits=8, $ xsize=xs*change, ysize=ys*change, $ xoff=(8.5-ys*change)/2.0, yoff=xs,$ /inches, /times, preview=1, encapsulated=1 endif else begin ; portrait device, file=postps, /portrait, bits=8, xsize=xs*changep, $ ysize=ys*changep, yoff=(11.0-ys*changep)/2.0, /color, $ xoff=(8.5-xs*changep)/2.0, /inches, /times, preview=1, encapsulated=1 endelse endelse endelse !P.THICK = 2.0 ; Slightly Thicker Line Weight for data lines assuming printer uses 300 dpi !X.THICK = 2.0 ; Slightly Thicker Line Weight for x-axis lines and tick marks !Y.THICK = 2.0 ; Slightly Thicker Line Weight for y-axis lines and tich marks endif else begin ; HP Laserjet print file set_plot, 'PCL', /copy, /interpolate ; color or no color if ycol eq 0.0 then begin ; no color if portrait_or_land_option eq 0 then begin ; landscape device, file=postps, /landscape, bits=8, xsize=xs*change, $ ysize=ys*change, yoff=11.0-(11.0-xs*change)/2.0, $ xoff=(8.5-ys*change)/2.0, /inches, /times endif else begin ; portrait device, file=postps, /portrait, bits=8, xsize=xs*changep, $ ysize=ys*changep, yoff=(11.0-ys*changep)/2.0, $ xoff=(8.5-xs*changep)/2.0,/inches, /times endelse endif else begin ; color if portrait_or_land_option eq 0 then begin ; landscape device, file=postps, /landscape, /color, bits=8, $ xsize=xs*change, ysize=ys*change, $ yoff=11.0-(11.0-xs*change)/2.0, xoff=(8.5-ys*change)/2.0, $ /inches, /times endif else begin ; portrait device, file=postps, /portrait, bits=8, xsize=xs*changep, $ ysize=ys*changep, yoff=(11.0-ys*changep)/2.0, /color, $ xoff=(8.5-xs*changep)/2.0, /inches, /times endelse endelse endelse ; endif else if abs(output_type) eq 4 then begin set_plot,'Z' ; Z-buffer device,set_resolution=[750,650] ; if plotting to the screen, the size is determined by the time length ; of the plot. This is because of the pixel size in the images. ; Postscript has adjustable pixel sizes, but the windows don't, so xsize ; must be adjusted to fit each time length endif else if output_type eq 3 then begin set_plot, 'X' device,pseudo_color=8 ; This IDL program is setup for 8bit graphics ; Note: window size is fixed window,1,xsize=750.,ysize=650.,colors=202 !P.FONT=-1 ; vector drawn font endif else if output_type eq -3 then begin ; MS Windows is IDL output device set_plot,'WIN' device,decomposed=0 ; Use 8-bit graphics on a 16-bit or 24-bit graphics system window,1,xsize=750.,ysize=650.,colors=202 ; Note: window size is fixed !P.FONT=-1 ; vector drawn font endif set_color_palette, output_type, ycol ; set up y tick marks for the electron and ion ranges: yticklab = ['10!A2','10!A3','10!A4'] ; Plot electron and ion energy spectrum ; if (strmid(partfile,0,4) eq 'null') or $ (strmid(partfile,0,4) eq 'NULL') $ then goto, cont6 if plotparsw ne 0 then print,' No problems found. Spectogram should be satisfactory.' plot_electron, ele2d, bt, et, nticks, tickval, blank, yticklab, $ change, plotparsw, dm_or_ssm_option,ycol plot_ion, ion2d, bt, et, nticks, tickval, blank, yticklab, $ change, flipi, plotparsw, dm_or_ssm_option,ycol if dm_or_ssm_option ne 4 and dm_or_ssm_option ne 5 then begin ; ; plot the electron and ion AVERAGE energy ; IN THE LOWER OF THE TWO B/W PLOTS ; ; plot the electron and ion INTEGRATED ENERGY FLUX ; or the INTEGRATED NUMBER flux according to the switch ; in the UPPERMOST OF THE TWO B/W PLOTS if plotparsw ne 0 then begin ; particle data found in interval if (intsw eq 1) then begin ; integrated energy flux plot_ee, eeflx, eave, ieflx, iave, blank, intsw, $ nticks, tickval, bt, et, yfil, yticklab, change endif else if (intsw eq 0) then begin ; integrated number flux plot_ee, oenflx, eave, oinflx, iave, blank, intsw, $ nticks, tickval, bt, et, yfil, yticklab, change endif endif endif ; cont6: ; Set up driftmeter tick marks if dm_or_ssm_option ne 3.0 then begin if dm_or_ssm_option eq 1 or dm_or_ssm_option eq 4 or $ dm_or_ssm_option eq 5 then begin if maxvel eq 0.0 then maxvel = 1.0 yb=-1.0*maxvel ye=maxvel endif else begin if max_nt eq 0.0 then max_nt = 1.0 yb=-1.0*max_nt ye=max_nt endelse ytickval=[yb,0.0,ye] yticklab=strcompress(string(fix(ytickval)),/remove_all) endif ; if there is no driftmeter file, we still have to plot the ephemeris ; data, but we use the J4/J5 file version, which has the mlt in seconds if (nodm eq 1) or (nodmdata eq 1) then begin mlt = mlt2 mlat = mlat2 endif ; if dm_or_ssm_option eq 1 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin print,'Begin plot of driftmeter data' plot_dm, dmdata, bt, et, nticks, tickval, ticklab, $ yb, ye, ytick, ytickval, yticklab, mlt, $ mlat, begtime, npts, cor, imul, nodm, ycol, $ change, vopt endif if dm_or_ssm_option eq 2 or dm_or_ssm_option eq 4 or dm_or_ssm_option eq 5 then begin ; We are using the J4/J5 values of mlat & mlt for consistency. ; These values are not found in the SSM data file but mlat(sub-satellite) ; is computed to find start and stop time and mlt could be computed. ; mlt = mlt2 mlat = mlat2 yb = -max_nt ye = max_nt ytickval=[yb,(yb+ye)/2.,ye] yticklab=strcompress(string(fix(ytickval)),/remove_all) ; print,'Begin plot of magnetometer data' plot_ssm,ssmdata, bt, et, nticks, tickval, $ yb, ye, ytick, ytickval, yticklab, $ begtime, nossm, ycol, $ change, output_type, dm_or_ssm_option endif if dm_or_ssm_option eq 4 then begin print,'Begin plot of SSIES/SM ion density data' plot_smni,smnidata, bt, et, nticks, tickval, ticklab,$ begtime, nosmni, $ change, ycol endif ; ; ; put the color scale on and the date and extra stuff plot_color_bar, cyb, cye, ctick, ctickval, cticke, cticki, $ change, diswitch, output_type, dm_or_ssm_option,ycol plot_extra, begtime, monname, satel, change, psfile, dm_or_ssm_option, $ bt, et, nticks, tickval, ticklab, imul, mlt, mlat, $ yb, ye, intsw, output_type, removebackground, JGF_LABEL ; ; go to the add_points procedure to see if there are lines or points to ; be added to the plot if n_elements(vcount) gt 0 then begin baset = float(begtime(3))*3600.0 + float(begtime(4))*60.0 add_points, vlines, pts, bt, et, vcount, pcount,$ vstr, vlabx, vlabc, change, baset, datafile_or_keyboard_input, $ dm_or_ssm_option endif ; we have completed the plotting, so close the file and set the plot ; type back to X-windows and exit the printing loop if (n_elements(outfile) gt 0) and $ (pagenum eq 1) and $ (first_time eq 1) then begin close, 50 openw, 50, outfile printf,50,'1 - Geographic Latitude in Degrees' printf,50,'2 - Geographic Longitude in Degrees' printf,50,'3 - Altitude in Km' printf,50,'4 - Invariant Latitude in Degrees' printf,50,'5 - Magnetic Local Time in Hours' if cor eq 1 then begin printf,50,'6 - Horizontal Velocity in Km/s' endif else begin printf,50,'6 - Horizontal Velocity - Corotation in Km/s' endelse printf,50,'7 - Vertical Velocity in Km/s' printf,50,'8 - Standard Deviation of Horizontal Velocity in Km/s' printf,50,'9 - Standard Deviation of Vertical Velocity in Km/s' printf,50,'10 - Incedent electron integral energy flux in '+ $ 'eV/(cm^2 s sr)' printf,50,'11 - Average incedent electron energy in eV' printf,50,'12 - Incedent electron integral number flux in '+ $ 'particles/(cm^2 s sr)' printf,50,'13 - Incedent ion integral energy flux in '+ $ 'eV/(cm^2 s sr)' printf,50,'14 - Average incedent ion energy in eV' printf,50,'15 - Incedent ion integral number flux in '+ $ 'particles/(cm^2 s sr)' printf,50,'16 - Electron integral energy flux at satellite'+ $ ' in eV/(cm^2 s sr)' printf,50,'17 - Average electron energy at satellite in eV' printf,50,'18 - Electron integral number flux at satellite'+ $ ' in particles/(cm^2 s sr)' printf,50,'19 - Ion integral energy flux at satellite'+ $ ' in eV/(cm^2 s sr)' printf,50,'20 - Average ion energy at satellite in eV' printf,50,'21 - Ion integral number flux at satellite'+ $ ' in particles/(cm^2 s sr)' printf,50,'SSYYMMDDHHMMSS 1 2 3 4'+ $ ' 5'+ $ ' 6 7 8 9'+ $ ' 10 11'+ $ ' 12 13 14 15'+ $ ' 16 17'+ $ ' 18 19 20 21' endif ; ; allow a break between x-window plots. ; ; ank = '' if abs(output_type) eq 3 then begin print,' Plot is over, to continue press any key.' Read, ank endif if abs(output_type) le 2 then begin ; PS mode print,'Spectogram plot is done' device, /close if output_type gt 0 then set_plot,'X' ; Switch from PS mode to X-Windows mode loop = 0 for i=0,strlen(postps)-1 do if strmid(postps,i,1) eq '.' then j=i filemain2 = strmid(postps,0,j) filext2 = strmid(postps,j,strlen(postps)-j) polarfile=filemain2+'_dm'+filext2 !p.font = 0 ; Graphics device specific font to be used if abs(output_type) eq 1 then set_plot, 'ps' else set_plot, 'PCL' ; device, file=polarfile, /landscape, /times if nodmdata eq 1 then nodm = 1 if nodmdata eq 1 then nodm = 0 device, /close if output_type gt 0 then set_plot,'x' endif else begin ; If we plotted to the screen, see if the user wants a hard copy if nodmdata eq 1 then nodm = 0 ; Transfer information from Z-buffer to file if abs(output_type) eq 4 then begin red = bytarr(256) green = bytarr(256) blue = bytarr(256) TVLCT, Red, Green, Blue, /GET outimage = bytarr(487500) ; create 750 x 650 pixel array ; get the image outimage = tvrd(0,0,750,650,0) year = begtime[0] month = begtime[1] dayofmonth = begtime[2] hour = begtime[3] minute = begtime[4] pngname = string(format='(a2,3i2.2,i4.4,2i2.2,a4)','cf',satel,$ month,dayofmonth, year, hour, minute,'.png') write_png,pngname,outimage, red, green, blue endif endelse first_time = 0 endif else print, 'J4/J5 Data does not exist for this interval.' endif ; increment the start time by the length of the plot loop_it: ;label used to skip plotting (liamk) begtime(4) = begtime(4) + length if begtime(4) ge 60 then begin begtime(3) = begtime(3) + 1 begtime(4) = begtime(4) - 60 endif ; calculate the "current time" in seconds curtime = begtime(3)*3600.0 + begtime(4)*60.0 endwhile if n_elements(outfile) gt 0 then close, 50 ; after we are all done with all of the times, close the j4 file j4cl = 1.0 ; THE END final1: end