;+ ; NAME: ; getipssources ; PURPOSE: ; Extract IPS observations from yearly data files in Morelia format ; CATEGORY: ; I/O ; CALLING SEQUENCE: FUNCTION getipssources, trange, $ observatory = observatory , $ subdir = subdir , $ use_good_sources = use_good_sources , $ degrees = degrees , $ silent = silent , $ point_sources = point_sources ; INPUTS: ; trange scalar or array[2]; type: time structure ; time range (begin and end time) ; a scalar is interpreted as [trange, trange+1 day] ; (i.e. a one day time range starting at trange ; OPTIONAL INPUT PARAMETERS: ; observatory scalar; type: string; default: none ; select one of 'mexart','nagoya', or 'ooty' ; (the list of valid observatories is hardcoded) ; subdir=subdir scalar; type: string; default: none ; subdirectory in $DAT/ips where to look for data files ; /degrees if set pp will be in degrees (default: radians) ; OUTPUTS: ; Result scalar; type: long; ; # ips sources in time range ; OPTIONAL OUTPUT PARAMETERS: ; point_sources=point_sources ; array[n]; type: structure ; structure with IPS sources inside time range ; point_sources.name = string with source name ; point_sources.pp = 2-element array with RA and decl of source ; point_sources.tt = time structure with time of observation ; point_sources.vv = IPS velocity ; point_sources.gg = IPS g-level ; point_sources.xy = used by PlotEarthSkymap to set pixel locations of sources ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; TimeFixYear, txt_read, TimeOp, TimeSet, TimeGet, TimeUnit ; CvSky, NewcombSun, elongation, BadValue ; CvPrecess, IsType, InitVar, vu_point_source ; PROCEDURE: ; > All relevant yearly IPS data files are read to find observations inside ; specified time range. These files are stored in $DAT/ips// ; where is supplied as keyword. ; > Files containing a list of IPS sources with celestial coordinates ; (presumably in B1950) are stored in directory $NAGOYA/sources ; A slightly different list is used each year. Currently this list ; contains the source name, B1950 right ascension and declination, ; and an expected source intensity. ; > RA and dec for the sources are precessed to the current epoch. ; MODIFICATION HISTORY: ; AUG-2016, Paul Hick (UCSD/CASS), using getnagoyasources as start ;- InitVar, use_good_sources, '' InitVar, silent , 0 observatories = ['mexart','nagoya','ooty'] IF IsType(observatory,/undefined) THEN message, 'no IPS observatory specified; select from '+strjoin(observatories,',') IF (where(observatory EQ observatories))[0] EQ -1 THEN message, 'invalid observatory specified; select from '+strjoin(observatories,',') dpm = ToDegrees(degrees=degrees) CASE n_elements(trange) OF 0: message, 'no time range specified in trange' 1: trange = [trange, TimeOp(/add,trange,TimeSet(/diff, day=1))] 2: IF IsType(trange,/string) THEN trange = TimeSet(trange) ELSE: trange = [trange[0], trange[n_elements(trange)-1]] ENDCASE jd_range = TimeGet(trange, /njd) yr_range = TimeGet(trange, TimeUnit(/year)) IF use_good_sources NE '' THEN yr_range[1] = yr_range[0] ;fmt = '(A8,I7,F6.2,F5.2,I4,I5,I4,I5,I5,I5,I4,F8.5)' ;fmt = '(A8,3I3,I4,2I3)' CASE IsType(subdir,/defined) OF 0: BEGIN InitVar, subdir, 'daily' ipsdir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='nagoya',subdir) END 1: BEGIN CASE CheckDir(subdir,/silent) OF 0: BEGIN InitVar, subdir, 'daily' ipsdir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='nagoya',subdir) END 1: ipsdir = subdir ENDCASE END ENDCASE srcdir = filepath(root=getenv('SSW_SMEI_DAT'),subdir='nagoya','sources') destroyvar, point_sources FOR yr=yr_range[0],yr_range[1] DO BEGIN ; Loop over yearly data files file = use_good_sources NE '' ? 'nagoya.'+use_good_sources+'.good' : 'nagoya.'+strcompress(yr,/rem) file = filepath(root=ipsdir,file) ; Construct file name ; Newer IPS data files have record length of 76 chars. The older ones ; sometimes do not have the scintillation index in the last column ; and have recordlength of 68 chars. status = txt_read(file, data, silent=silent, /test) ; Read yearly IPS file IF NOT status THEN status = txt_read(file, data, approx_recl=68, silent=silent) ; Read yearly IPS file IF status THEN BEGIN srcs = strmid(data,0,8) yrs = TimeFixYear( fix( strmid(data,9,2) ) ) mons = fix ( strmid(data,11,2) ) days = fix ( strmid(data,13,2) ) uts = float( strmid(data,15,6) ) vvs = float( strmid(data,49,5) ) ;ras = AngleUnits( from_time=strmid(data,69,9), /to_degrees) ;decs = AngleUnits( from_sexa=strmid(data,79,9), /to_degrees, /singlesign ) ; files with ACE data are mangled tqus = strmid(data,54,4) iqus = where(strlen(strcompress(tqus,/remove)) eq 0) IF iqus[0] NE -1 THEN tqus[iqus] = '0' qus = fix(tqus) ;strlen(strcompress(tqus,/remove)) eq 0) ? 0 : fix(tqus) ; qus = fix ( strmid(data,54,4) ) ggs = float( strmid(data,58,8) ) ; Check for 'bad data'. For velocities check for value -999 or ; quality indicator -999. For g-level check for value 0.0 tmp = where(vvs EQ -999.0) ; OR qus EQ -999) IF tmp[0] NE -1 THEN vvs[tmp] = BadValue(vvs) tmp = where(ggs EQ 0.0) IF tmp[0] NE -1 THEN ggs[tmp] = BadValue(ggs) tts = TimeSet(yr=yrs, mon=mons, day=days, hour=uts) ; Time of obs jds = TimeGet(tts, /njd) ; Time of obs in Julian days ;isrc = where(jd_range[0] LE jds AND jds LE jd_range[1] AND qus NE -999, ni) ; Obs inside time range jd_range ; Pick up all sources inside time range with either V or g good isrc = where(jd_range[0] LE jds AND jds LE jd_range[1] AND (finite(vvs) OR finite(ggs)), ni) IF ni NE 0 THEN BEGIN ; If obs inside time range found pps = replicate(BadValue(0.0),2,ni) file = filepath(root=srcdir,'nagoyasources.'+strcompress(yr,/rem)) ; List of obs sources for relevant year ; recl spec should not be necessary here (just leads ; to a call to the crappy GuessRecs). ;IF txt_read(file, data, approx_recl=39, silent=silent) THEN BEGIN ; Read file with source positions IF txt_read(file, data, silent=silent) THEN BEGIN ; Read file with source positions ips_src = strmid(data,0,8) ra = AngleUnits( from_time=strmid(data, 8,9), /to_degrees) dec = AngleUnits( from_sexa=strmid(data,18,9), /to_degrees, /singlesign ) FOR i=0,ni-1 DO BEGIN j = where(srcs[isrc[i]] EQ ips_src) CASE j[0] NE -1 OF 0: say, 'no RA,dec for IPS source '+srcs[isrc[i]], threshold=1, silent=silent 1: pps[*,i] = [ra[j[0]],dec[j[0]]] ENDCASE ENDFOR ENDIF i = vu_point_source(i, srcs[isrc], /name, create=ni) ; Create the vu_point_source structure i i = vu_point_source(i, pps , /location, /degrees) ; Update structure i = vu_point_source(i, tts[isrc] , /time) i = vu_point_source(i, vvs[isrc] , /ips_speed ) i = vu_point_source(i, ggs[isrc] , /ips_glevel ) boost, point_sources, i ENDIF ENDIF ENDFOR ; Precess from B1950 to current epoch IF IsType(point_sources, /struct) THEN BEGIN tmp = CvPrecess(TimeSet(bepoch=1950.0d0), $ vu_point_source(point_sources,/time), $ vu_point_source(point_sources,/location,degrees=degrees), $ degrees=degrees) point_sources = vu_point_source(point_sources, tmp, /location, degrees=degrees) ENDIF dT = TimeOp(/subtract,trange[1],trange[0],TimeUnit(/hour)) IF dT/2*2 EQ dT THEN dT /= 2 ELSE dT /= 2.0 TT = TimeOp(/meant,trange[0],trange[1]) nsources = n_elements(point_sources) say, strcompress(nsources) +' IPS sources in '+ $ TimeGet(TT,/ymd,upto=TimeUnit(/minute))+' [+/-'+strcompress(dT)+' hr]', threshold=1, silent=silent RETURN, nsources & END