;+ ; NAME: ; nagoya_glevel ; PURPOSE: ; Estimate G-level from scintillation index ; CATEGORY: ; sat/idl/toolbox/ips ; CALLING SEQUENCE: PRO nagoya_glevel, TT, $ plotx = plotx , $ all_data = all_data , $ input = input , $ output = output , $ rnum = rnum , $ nsg = nsg , $ whole_year = whole_year , $ use_min_glevel = use_min_glevel, $ min_ppdist = min_ppdist , $ _extra = _extra ; INPUTS: ; TT array[1]; type: time structure ; only g-levels for source observation in TT +/- 0.5 days ; are updated (unless /whole_year is set). ; OPTIONAL INPUT PARAMETERS: ; input=input scalar; type: string; default: $SSW_SMEI_DAT/nagoya/daily ; fully-qualified name of a yearly Nagoya IPS file of type ; nagoya.yyyy, or directory where these files are located. ; If a directory is specified then the most recent yearly file ; is used. ; The default is the directory that holds the yearly files ; derived from the real-time Nagoya IPS data. These are the ; files created and maintained by dailyips.f by ingesting ; the VLIST_UCSD_* files downloaded from Nagoya daily. ; /all_data if set, g-level is calculated from all data. ; if not, g-level is calculated from past data. ; /use_min_glevel if scintillation indices are available from multiple stations ; calculated g-level for all stations and use the lowest ; g-level as the final value. ; ; /output if set, update the input file ; output=output scalar; type: string; the name of a file into which the ; updated data are written (instead of overwriting the input ; file if /output is specified). ; ; ; /plotx if set, a GIF map is created and displayed ; rnum=rnum scalar; type: integer; default: 8 ; minimum # data point required for reliable fitting ; nsg=nsg scalar; type: float; default: 1.5 ; criterion to eliminates highly-scattered data ; from fitted line. (nsg*sigma) ; /whole_year if set, all g-levels in the yearly file are recalculated ; (argument TT, if set, is ignored). ; min_ppdist=min_ppdist ; scalar; type: float; default: 0.2 ; minimum point-P distance in AU. IPS observations with distances ; smaller than this are not used in the calculation of the scint index ; vs. point-P distance dependence of each IPS source. ; OUTPUTS: ; If /plotx is set then the sources in TT +/- 0.5 days are displayed in a fisheye plot. ; If keyword 'output' is used then the modified G-levels are written to file. ; If neither keyword is used then a 'test run' is made with only text messages ; to the screen. ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; IsType, txt_read, TimeGet, TimeFixYear, IsTime, TimeUnit ; TimeGet, TimeSet, GetColors, CvPointOnLos, FishEye, twin, boost, InitVar ; flt_string ; PROCEDURE: ; The IPS data should be stored in yearly files with ; name 'nagoya.yyyy' where yyyy is the year. ; ; IPS DATA FORMAT: (OLD) ; 0 1 2 3 4 5 6 7 ; 0123456789012345678901234567890123456789012345678901234567890123456789012345 ; SOURCE YRMNDY UT DIST HLA HLO GLA GLO CARR V ER G-LEVEL SC-INDX G-LEVEL SC-INDX G-LEVEL SC-INDX ; 1117+14 000711 7.10 0.85 6 -33 6 316 1965 303-999 0.00000 0.518E+02 0.00000 0.518E+02 0.00000 0.518E+02 ; ; The second pair g-level/scint index pair contains the g-level and scint index for Kiso. ; The third pair contains the g-level and scint index for Fuji. ; ; The first g-level/scint index pair is a 'best value' pair. Currently the Fortran program ; href=dailyips= puts the scint index for Kiso there, while the g-level gets filled ; in here depending on keyword selection. By default the Kiso g-level is used; ; if /use_min_glevel is set then the smallest g-level (Kiso or Fuji) is used. ; ; Scintilation index values, S, are fitted to a power law function of the point-P distance, r: ; S = 10^a * r^b ; (a,b) are determined by a linear least-square fit. ; The g-level, g, for an individual scintillation index follows by normalizing to the ; fitted function: ; g = sqrt( S/(10^a+r^b)) ; Note that the g-level for a fixed scintillation will change if the fit changes (i.e. if ; a different set of observations is used to determine the fit. ; MODIFICATION HISTORY: ; JUL-2000, K. Fujiki (STELAB) ; JUL-2000, Paul Hick (UCSD/CASS) ; The calculation of g-levels was modified in two ways: ; 1. The calculation of standard deviations is now done with the IDL ; routine 'stddev', rather than Fujiki's function 'sigma' (which gave ; a slightly different value) ; 2. The sources for which the g-levels are to be recalculated are all sources ; within 0.5 days of the specified time. Fujiki's original code would sometimes ; miss a few sources when the time would be close to the start or end of a month. ; The 'fish-eye' plot now shows the elongation in the radial direction (consistent ; with the display in href=vu_earthskymap= (Fujiki's original code used the point-P ; distance as the radial coordinate). ; SEP-2001, Paul Hick (UCSD/CASS) ; Added the /whole_year keyword to process all observations in a yearly file in one pass. ; SEP-2002, Paul Hick (UCSD/CASS) ; Modified to handle files with multiple scintilation indices. ; Two modifications were required: the txt_read call now uses ; /test to determine the record length instead of a hardcoded length ; of 76. The output section truncated records after the first ; scintillation index. Now it doesn't. The current modifications ; should allow additions of other scint indices in the future. ; NOV-2002, Paul Hick (UCSD/CASS) ; Added separate calculations of Kiso and Fuji g-levels. ; Added keyword /use_min_glevel (which determines calculation of a 'best' g-level. ; JUL-2003, Paul Hick (UCSD/CASS) ; Fixed bug in fitting procedure for g-levels (in nagoya_r_fitting). It would crash ; if only two data values are available. ; APR-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Fixed bug in TimeSet with /botime set (keywords changed some time ago) ; AUG-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Added some documentation, Put two helper functions into separate files. ;- InitVar, all_data , /key InitVar, plotx , /key InitVar, whole_year , /key InitVar, use_min_glevel , /key InitVar, rnum , 8 ; Reliability criterion for eliminating InitVar, nsg , 1.5 ; .. scattered data point InitVar, min_ppdist , 0.2 ; Point-P distances are stored in F4.2 format. ; To avoid problems with log(0.0) in the fit between log(scint indx) and log(point-P dist) ; set the min_ppdist at least to 0.009. min_ppdist >= 0.009 use_best_station = 1-use_min_glevel ; The default directory is $SSW_SMEI_DAT/nagoya/daily. Input files are named nagoya.yyyy InitVar, input, filepath(root=getenv('SSW_SMEI_DAT'), subdir=['nagoya','daily'], '') CASE checkdir(input,/silent) OF 0: BEGIN files = file_search(input, count=icnt) IF icnt EQ 0 THEN BEGIN message, /info, 'file not found: '+input RETURN ENDIF END 1: BEGIN file19 = filepath(root=input, 'nagoya.19*') file20 = filepath(root=input, 'nagoya.20*') files19 = file_search(file19, count=icnt19) files20 = file_search(file20, count=icnt20) icnt = icnt19+icnt20 IF icnt EQ 0 THEN BEGIN message, /info, 'file not found: '+file19+' '+file20 RETURN ENDIF IF icnt19 GT 0 THEN boost, files, files19[0:icnt19-1] IF icnt20 GT 0 THEN boost, files, files20[0:icnt20-1] END ENDCASE ; Work with the most recent yearly file only input_file = files[icnt-1] ; Check for output. If keyword output is set to a string it is assumed to ; the name of the output file. If not then the input file (most recent ; yearly file is overwritten). CASE IsType(output, /string) OF 0: BEGIN output_file = input_file InitVar, output, /key END 1: BEGIN output_file = output output = 1 END ENDCASE ; Read data from yearly file as string array. IF NOT txt_read(input_file, cobs, nrec=nobs, /test) THEN $ RETURN ut = TimeSet(yr=TimeFixYear(long(strmid(cobs,9,2) )), $ month = long ( strmid(cobs,11,2) ), $ day = long ( strmid(cobs,13,2) ), $ hour = float( strmid(cobs,16,5) ) ) jd = TimeGet(ut, /njd) ; If no time is specified, take the last time on file ut0 = IsTime(TT) ? TT[0] : TimeGet( /botime, ut[nobs-1], TimeUnit(/day) ) jd0 = (TimeGet(ut0, /njd))[0] source = strtrim(strmid(cobs, 0,9)) ; Source name ppdist = float (strmid(cobs,22,4)) ; Point P distance (AU) ; G-level/scintillation index pairs start at column 58: one 'best pair', followed by ; pairs for individual stations, currently (2002) Kiso and Fuji. ; Each pair consists of 8+10=18 chars in (F8.5,E10.3) format. There should always ; be a blank space at columns 58 and 66 (the first char of each number). glevel_start = 58 ; !!!!!!!! FUDGE FUDGE FUDGE !!!!!!!!!!! ; ; Keep g-level between glevel_min and glevel_max. ; This fits the format F8.5 used for the g-levels while leaving a blank space at the start. ; A g-level outside this range is set to glevel_bad glevel_min = 0.0 ; Minimum allowed g-level glevel_max = 9.99 ; Maximum allowed g-level glevel_bad = 0.0 ; No g-level available scindx = strmid(cobs[0],glevel_start) ; Columns with g-level/scintillation index pairs scindx = flt_string(scindx, /exp, fmt=fmt) ; Get format from first record nstation = n_elements(scindx)/2 ; # scintillation indices per record use_best_station OR= nstation LE 2 station_name = ['best','Kiso','Fuji'] best_station = 0 CASE use_best_station OF 0: BEGIN first_station = 1 last_station = nstation-1 END 1: BEGIN first_station = best_station last_station = best_station END ENDCASE ipsdata = fltarr(2*nstation, nobs) ; Get all scint index values reads, strmid(cobs,glevel_start), format=fmt, ipsdata glevel = 0+2*indgen(nstation) ; g-level postions in ipsdata scindx = glevel+1 ; scint index positions in ipsdata ; Process either the best index only, or loop over all available stations. FOR station=first_station,last_station DO BEGIN glevel_val = reform(ipsdata[glevel[station],*]) scindx_val = reform(ipsdata[scindx[station],*]) ; Select data for specified day (usually the latest day) ; Get list of indices 'inew' into 'cobs' array. These are the sources for ; which the g-level is updated. CASE whole_year OF 0: BEGIN djd = 0.5 inew = where(jd0-djd LT jd AND jd LT jd0+djd AND scindx_val GT 0, icnt) END 1: BEGIN inew = uniq(source,sort(source)) icnt = n_elements(inew) message, /info, 'total # sources:'+strcompress(icnt) END ENDCASE CASE icnt OF 0: message, /info, station_name[station]+' station: no sources in '+TimeGet(ut0, /ymd)+' +/- 12 hours' ELSE: BEGIN ; Calculate g-levels FOR i=0L,icnt-1 DO BEGIN ; Loop over sources ; Indices to valid observations for source i ; The test for positive ppdist is not really necessary. Sometime ppdist=0.0 for observations ; very close to the Sun, causing problems when log10(ppdist) is calculated. imod = where(source EQ source[inew[i]] AND ppdist GT 0.0 AND scindx_val GT 0.0, inum) isrc = imod IF inum GE rnum AND NOT all_data THEN BEGIN ; Real-time update: use only past observations a = where(jd[isrc] LT jd0+0.5, inum) IF inum GE rnum THEN isrc = isrc[a] ENDIF IF inum GE rnum THEN BEGIN ; Enough data for fitting a = nagoya_r_fitting(ppdist[isrc],scindx_val[isrc],min_ppdist) IF a[1] LT 0 THEN BEGIN ; Good fit ; Difference in LOG space dif = alog10(scindx_val[isrc]) - (a[1]*alog10(ppdist[isrc])+a[0]) ; Eliminate outliers a = where(abs(dif) LT nsg*stddev(dif), inum) IF inum GE rnum THEN BEGIN isrc = isrc[a] ; Fit again after removing outliers a = nagoya_r_fitting(ppdist[isrc],scindx_val[isrc],min_ppdist) glevel_val[imod] = sqrt(scindx_val[imod]/(10^(a[1]*alog10(ppdist[imod])+a[0]))) ENDIF ENDIF ENDIF ENDFOR IF NOT whole_year THEN $ message, /info, station_name[station]+' station: '+strcompress(icnt, /rem) +' IPS sources in '+ $ TimeGet(ut0,/ymd,upto=TimeUnit(/minute))+' [+/-'+strcompress(round(djd*24))+' hr]' a = where(ipsdata[glevel[station],*] NE glevel_val, icnt) message, /info, station_name[station]+' station: '+strcompress(icnt,/rem)+'/'+strcompress(nobs,/rem)+' g-levels modified' ; Flag g-level outside range [glevel_min, glevel_max] as bad a = where(glevel_val LT glevel_min OR glevel_val GT glevel_max) IF a[0] NE -1 THEN BEGIN message, /info, 'discarded negative and very large g-levels' glevel_val[a] = glevel_bad ENDIF ipsdata[glevel[station],*] = glevel_val END ENDCASE ENDFOR IF use_min_glevel THEN BEGIN ; Set g-level to minimum value from all real stations glevel_val = reform(ipsdata[glevel[best_station],*]) FOR i=0L,nobs-1 DO BEGIN tmp = ipsdata[glevel[first_station:last_station],i] a = where(tmp GT 0) IF a[0] NE -1 THEN tmp = min(tmp[a]) ELSE tmp = 0.0 ipsdata[glevel[best_station],i] = tmp ENDFOR a = where(glevel_val ne ipsdata[glevel[best_station],*], icnt) message, /info, station_name[best_station]+' station: '+strcompress(icnt,/rem)+'/'+strcompress(nobs,/rem)+' g-levels modified' ENDIF IF output THEN BEGIN ; Output section ; The loop is used because the string function truncates at 1024 lines. a = strarr(nobs) FOR n=0L,nobs-1,1024 DO a[n: (n+1023) < (nobs-1)] = string(ipsdata[*,n: (n+1023) < (nobs-1) ],format=fmt) a = strmid(cobs,0,glevel_start)+a openw, /get_lun, iu, output_file printf, iu, a, format='(A)' free_lun, iu message, /info, 'new G-levels to '+output_file ENDIF IF plotx THEN BEGIN ; Plotting section glevel = reform(glevel[best_station,*]) rp = float ( [ [(strmid(cobs,31,4))[inew]], $ ; Longitude [(strmid(cobs,26,4))[inew]] ] ) ; Latitude rp = transpose( [ [rp/!radeg], [ppdist[inew]] ] ) ; Add point-P distance vel = float( (strmid(cobs,50,4))[inew] ) ; Velocity twin, /show, xs=500,ys=500 !p.color = !d.n_colors-1 !p.background = 0 stretch, !d.n_colors-1, 0 erase xyouts, 0.5, 0.95, TimeGet(ut0,/ymd), /normal, si=1.5, align=0.5 nagoya_plotg2d, rp, glevel[inew], source=source[inew], $ color=GetColors(vel, 100+100*indgen(10), /legend, /open, _extra=_extra) ; flip_colors ENDIF RETURN & END