;+ ; NAME: ; vu_header ; PURPOSE: ; Reads header from nv3d* or wson* file ; CATEGORY: ; I/O ; CALLING SEQUENCE: FUNCTION vu_header, files, check=check, index=index, count=count, $ sort_time=sort_time, allow_nrad_changes=allow_nrad_changes, silent=silent ; INPUTS: ; files array; type: string ; list of file names ; OPTIONAL INPUT PARAMETERS: ; /check forces consistency checks ; - check array dimensions ; - check source surface distance and radial separation ; - check data type ; /sort_time sort headers into chronological order before returning ; (if the 'files' array is output from vu_select this is done ; already) ; /allow_nrad_changes ; if /check is set this bypasses the check on the number of ; radial distances. The start distance and step size must be ; identical for all headers, but the number of distances ; can vary (href=vu_read= will only read the minimum number ; of distances specified in all headers). ; OUTPUTS: ; Result array; type: vu_header structure ; header information for all files ; if /sort_time is set then the headers are sorted ; into chronological order before returning. ; files array; type: string ; input array with files that are not ; valid nv3d* or wson* files remove ; OPTIONAL OUTPUT PARAMETERS: ; count scalar; type: long integer ; # entries returned in arrays 'files' and 'hdrs' ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, flt_string, BadValue, destroyvar, TimeSet, vu_check ; flt_read, SetFileSpec, GetFileSpec, TimeOp ; EXTERNAL: ; vu_header__define ; RESTRICTIONS: ; Using the /sort_time keyword may have unexpected effects when reading files ; with markers in the file name (from the time-dependent tomography) due to ; precision problems with the Carrington variables used in the file name ; (the last digit is just barely significant). ; PROCEDURE: ; Lines in the file header are used to fill a structure containing the ; information about the file content. The exact strings used to identify the ; information in each line should match the strings used in the Fortran ; subroutine href=vu_header=. ; MODIFICATION HISTORY: ; FEB-2001, Paul Hick (UCSD/CASS) ; JAN-2001, Paul Hick (UCSD/CASS) ; Added /exp to flt_string call for reading powers for density and distance ; dependence of density flucuations ; JUL-2002, Paul Hick (UCSD/CASS) ; V1.03; added nl, nlos, dlos and binx ; JUN-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Added capability to handle *.gz files transparently ;- InitVar, check , /key InitVar, sort_time , /key InitVar, allow_nrad_changes, /key InitVar, silent , 0 count = n_elements(files) hdrs = replicate({vu_header}, count) FOR i=0,count-1 DO BEGIN CASE flt_read(files[i], comments=hdr, /header,silent=2) OF ; Read header 0: IF silent LE 1 THEN message, /info, 'error accessing '+hide_env(files[i]) 1: BEGIN hdr = strmid(hdr,2) pos = (where(strpos(hdr,'created on') NE -1))[0] ; Contains file name and time of file creation IF pos NE -1 THEN BEGIN str = strtok(hdr[pos], /extract) SetFileSpec, files[i] ; Fully-qualified file name hdrs[i].file = GetFileSpec(upto=(['type','name'])[GetFileSpec(part='type') EQ '.gz']) t = n_elements(str) t = TimeSet(str[t-2]+' '+str[t-1]) hdrs[i].created = t ; Creation time ENDIF ; Pick up value used in file to indicate bad data points (these will be ; replaced by NaN numbers by vu_read) pos = (where(strpos(hdr,'Bad value flag:') NE -1))[0] IF pos NE -1 THEN hdrs[i].bad = flt_string(hdr[pos], /double) ELSE hdrs[i].bad = BadValue(hdrs[i].bad) ; The presence of a 'T3D_header, v' line (containing a version number) is used ; to validate the file as a valid tomography file. pos = (where(strpos(hdr,'T3D_header, v') NE -1))[0] IF pos EQ -1 THEN $ message, /info, 'not a tomography file: '+hide_env(files[i]) $ ELSE BEGIN ; Pick up version number of header. IF pos NE -1 THEN hdrs[i].version = strmid(hdr[pos], strlen('T3D_header, v')) ; The Carrington offset is stored as a separate number in the file header. ; Here we pick it up, but do not store in the header structure. Instead it ; gets added to all the appropriate Carrington variables. cst = 'File name prefix:' pos = (where(strpos(hdr,cst) ne -1))[0] CASE pos NE -1 OF 0: BEGIN SetFileSpec, hdrs[i].file hdrs[i].prefix = strmid(GetFileSpec(part='name'),0,4) END 1: hdrs[i].prefix = strcompress(strmid( hdr[pos], strlen(cst) ), /rem) ENDCASE pos = (where(strpos(hdr,'Carrington offset:') NE -1))[0] IF pos NE -1 THEN offset = flt_string(hdr[pos], /double) cst = 'Carrington time:' pos = (where(strpos(hdr,cst) NE -1))[0] IF pos NE -1 THEN hdrs[i].uttime = vu_set_time_entry(hdr[pos], cst, offset) pos = (where(strpos(hdr,'Carrington limits of array (start/end):') NE -1))[0] IF pos NE -1 THEN hdrs[i].array = offset+flt_string(hdr[pos], /double) pos = (where(strpos(hdr,'Carrington limits of region of interest (start/end):') NE -1))[0] IF pos NE -1 THEN hdrs[i].roi = offset+flt_string(hdr[pos], /double) pos = (where(strpos(hdr,'Carrington start time:') NE -1))[0] IF pos NE -1 THEN hdrs[i].time_start = offset+flt_string(hdr[pos], /double) pos = (where(strpos(hdr,'Carrington time resolution:') NE -1))[0] IF pos NE -1 THEN hdrs[i].time_step = flt_string(hdr[pos], /double, /exp) cst = 'Carrington forecast time:' pos = (where(strpos(hdr,cst) NE -1))[0] IF pos NE -1 THEN hdrs[i].time_forecast = vu_set_time_entry(hdr[pos], cst, offset) pos = (where(strpos(hdr,'Latitude range (degrees):') NE -1))[0] IF pos NE -1 THEN hdrs[i].latitude_range = flt_string(hdr[pos], /exp) ELSE hdrs[i].latitude_range = [-90.0,90.0] pos = (where(strpos(hdr,'Radial reference distance (AU):') NE -1))[0] IF pos NE -1 THEN hdrs[i].distance_start = flt_string(hdr[pos], /exp) pos = (where(strpos(hdr,'Radial resolution (AU):') NE -1))[0] IF pos NE -1 THEN hdrs[i].distance_step = flt_string(hdr[pos], /exp) pos = (where(strpos(hdr,'Power index of density dependence of density fluctuations (V-data,G-data):') NE -1))[0] IF pos NE -1 THEN hdrs[i].power_density = flt_string(hdr[pos], /exp) pos = (where(strpos(hdr,'Power index of radial dependence of density fluctuations (V-data,G-data):') NE -1))[0] IF pos NE -1 THEN hdrs[i].power_distance = flt_string(hdr[pos], /exp) pos = (where(strpos(hdr,'Density at 1 AU (cm^-3):') NE -1))[0] IF pos NE -1 THEN hdrs[i].density_at_1AU = (flt_string(hdr[pos]))[2] pos = (where(strpos(hdr,'Spatial filters for smoothing velocity and density (degrees):') NE -1))[0] IF pos NE -1 THEN hdrs[i].smooth_filter = flt_string(hdr[pos]) pos = (where(strpos(hdr,'Spatial filters for filling velocity and density (degrees):') NE -1))[0] IF pos NE -1 THEN hdrs[i].fill_filter = flt_string(hdr[pos]) pos = (where(strpos(hdr,'Temporal filters for velocity and density (days):') NE -1))[0] IF pos NE -1 THEN hdrs[i].time_filter = flt_string(hdr[pos]) pos = (where(strpos(hdr,'Clip longitude (degrees):') NE -1))[0] IF pos NE -1 THEN hdrs[i].clip_longitude = flt_string(hdr[pos]) pos = (where(strpos(hdr,'Dimensions (nLng,nLat,nRad,nTim):') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].nLng = tmp[0] hdrs[i].nLat = tmp[1] hdrs[i].nRad = tmp[2] hdrs[i].nTim = tmp[3] ENDIF pos = (where(strpos(hdr,'Iteration:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].iteration = tmp[0] hdrs[i].niterations = tmp[1] ENDIF pos = (where(strpos(hdr,'Time index:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].time_index = tmp[0] ENDIF pos = (where(strpos(hdr,'# Velocity lines of sight:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].nl[0] = tmp[0] ENDIF pos = (where(strpos(hdr,'# G-level lines of sight:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].nl[1] = tmp[0] ENDIF pos = (where(strpos(hdr,'# Line-of-sight segments for velocity:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].nlos[0] = tmp[0] ENDIF pos = (where(strpos(hdr,'# Line-of-sight segments for g-level:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].nlos[1] = tmp[0] ENDIF pos = (where(strpos(hdr,'Line-of-sight resolution for velocity:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].dlos[0] = tmp[0] ENDIF pos = (where(strpos(hdr,'Line-of-sight resolution for g-level:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].dlos[1] = tmp[0] ENDIF pos = (where(strpos(hdr,'# segments/bin for velocity:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].binx[0] = tmp[0] ENDIF pos = (where(strpos(hdr,'# segments/bin for density:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = round(flt_string(hdr[pos])) hdrs[i].binx[1] = tmp[0] ENDIF pos = (where(strpos(hdr,'Linear scaling constants:') NE -1))[0] IF pos NE -1 THEN BEGIN tmp = flt_string(hdr[pos], /exp) CASE n_elements(tmp) OF 2 : hdrs[i].scale = [tmp,tmp] ELSE: hdrs[i].scale = tmp ENDCASE ENDIF pos = (where(strpos(hdr,'Power for radial normalization:') NE -1))[0] IF pos NE -1 THEN hdrs[i].r_pwr = flt_string(hdr[pos]) pos = (where(strpos(hdr,'Rotation counter:') NE -1))[0] IF pos NE -1 THEN hdrs[i].marker = round(flt_string(hdr[pos])) pos = (where(strpos(hdr,'Velocity (km s^-1)/normalized density (cm^-3)') NE -1))[0] hdrs[i].plasma = pos NE -1 pos = (where(strpos(hdr,'Radial (r^2*B_r)/tangential (r*B_t) magnetic field') NE -1))[0] hdrs[i].magnetic = pos NE -1 pos = (where(strpos(hdr,'User:') NE -1))[0] if pos ne -1 then hdrs[i].user_info = strmid(hdr[pos], strlen('User:')+1) ENDELSE END ENDCASE ENDFOR ; Discard all files which are not recognized as valid t3d* or b3d* files index = indgen(count) tmp = where(hdrs.version NE '', count) IF count NE 0 THEN BEGIN hdrs = hdrs [tmp] files = files[tmp] index = index[tmp] ENDIF ELSE BEGIN destroyvar, files, index hdrs = {vu_header} ; Return blank header structure ENDELSE IF count NE 0 AND check THEN BEGIN ; Checks are done by comparing with the first header, so at least one header ; (the first one) will remain. nLng = hdrs[0].nLng nLat = hdrs[0].nLat nRad = hdrs[0].nRad ; Check dimensions. If /allow_nrad_changes is set don't check nrad. CASE allow_nrad_changes OF 0: tmp = hdrs.nLng EQ nLng AND hdrs.nLat EQ nLat AND hdrs.nRad EQ nRad 1: tmp = hdrs.nLng EQ nLng AND hdrs.nLat EQ nLat ENDCASE ntmp = round(total(tmp)) IF ntmp NE count THEN BEGIN FOR i=0,count-1 DO IF NOT tmp[i] THEN message, /info, 'dump '+hdrs[i].file+', bad dims' tmp = where(tmp) hdrs = hdrs [tmp] files = files[tmp] index = index[tmp] count = ntmp ENDIF ; Check source surface distance and radial separation tmp = hdrs.distance_start[0] EQ hdrs.distance_start AND $ hdrs.distance_step [0] EQ hdrs.distance_step ntmp = round(total(tmp)) IF ntmp NE count THEN BEGIN FOR i=0,count-1 DO IF NOT tmp[i] THEN message, /info, 'discarding '+hdrs[i].file+': not the radial grid' tmp = where(tmp) hdrs = hdrs [tmp] files = files[tmp] index = index[tmp] count = ntmp ENDIF tmp = hdrs.plasma [0] EQ hdrs.plasma AND $ hdrs.magnetic[0] EQ hdrs.magnetic ntmp = round(total(tmp)) IF ntmp NE count THEN BEGIN FOR i=0,count-1 DO IF NOT tmp[i] THEN message, /info, 'discarding '+hdrs[i].file+': wrong data type' tmp = where(tmp) hdrs = hdrs [tmp] files = files[tmp] index = index[tmp] count = ntmp ENDIF ENDIF IF count NE 0 AND sort_time THEN BEGIN tmp = sort( TimeOp(/subtract,hdrs.uttime, hdrs[0].uttime, TimeUnit(/day)) ) hdrs = hdrs [tmp] files = files[tmp] index = index[tmp] ENDIF RETURN, hdrs & END