; Comments added by Erik Johnson 7/16/07
; read_maccs.pro, compiles and runs with idl
; called by igloolik_despike.pro (which creates despiked .s2 files)
; called by five_sec_avg.pro (which creates 5 second average .5s files)
; called by sec_to_min.pro (which creates 1 minute average .1m files)
; reads 5 different types of maccs data: raw maccs, cleaned maccs, raw gsc, cleaned gsc, 5 sec & 1 min
; GSC = Geological Survey of Canada (The Geomagnetic Monitoring Service of the GSC operates a network of 13 Magnetic Observatories across Canada which includes all CANMOS magnetometers and the Alert magnetometer. It used to operate the Mould Bay magnetometer which is now closed.)
; calculates and returns time
; stores the data in variables and returns them
; the only required parameter is the filename unless it is a raw file (then raw must be set to 1 or 2)
; Ex: read_maccs, "CD04025.s2"
; if you set the diag parameter (diag="something") it outputs the following to the command line:
; *****************************************************************************
; if you are reading a raw file you must set the raw parameter to 1 for a maccs file or 2 for a gsc file.
; to prevent the program from printing the name of the file it is processing: quiet="something"
; If you set the nt parameter (ex: nt="nt") the field values will be divided by 1000 to get the values in nanoteslas.
; Ex: read_maccs, "/Volumes/physics_data/Ftp/MACCS_DATA/Raw/CD/2007ir/CD07193.dat", diag="bob", raw=1, nt="nanoT", quiet="please"
pro read_maccs,filename,time,bx1,by1,bz1,type,raw=raw,nt=nt,$
bx2=bx2,by2=by2,bz2=bz2,tf1=tf1,int=int,diag=diag,quiet=quiet,block=block
; set the keyword raw to 1 for raw maccs and to 2 for raw gsc
; x and y are the byte dimensions of the file. The order for the different file
; types is raw maccs, raw gsc, cleaned maccs, cleaned gsc, 5 sec, and 1 min
; (the raw files do not have fixed length)
x = [35, 21, 28, 15, 15, 14]
y = [00, 00, 86400L, 86400L, 17280L, 1440L]
sizes = [00, 00, 2419200L, 1296000L, 259200L, 20160L]
; Open file and get file size
if not keyword_set(quiet) then print,'Now processing ',filename
openr,unit,filename,error=err,/get_lun
if err ne 0 then print,'no file ',filename
file_info = fstat(unit)
file_size = file_info.size
if keyword_set(raw) then begin
type = raw - 1
if raw eq 1 then y(0) = file_size / x(0)
if raw eq 2 then y(1) = file_size / x(1)
endif else begin
type = where(file_size eq sizes, count)
if count eq 0 then begin
print,'Error in file size.'
stop
endif else type = type(0)
endelse
if (err eq 0) then begin
m=x(type)
n=y(type)
data_array = bytarr(m,n)
readu,unit,data_array
free_lun,unit
endif else begin
print,'An error has occurred while trying to read.'
stop
endelse
; Extract the individual data arrays from the file byte array. Also, divide
; the field values by 1000. to get nanotesla.
; all files
offset = [3, 2, 1, 0, 0, 0]
hh = data_array(offset(type) + 0,*)
mm = data_array(offset(type) + 1,*)
offset(5) = -1
bx1 = long(data_array(offset(type) + 3:offset(type) + 6,*),0,1,y(type))
by1 = long(data_array(offset(type) + 7:offset(type) + 10,*),0,1,y(type))
bz1 = long(data_array(offset(type) + 11:offset(type) + 14,*),0,1,y(type))
case type of
5: ss = 0
4: ss = data_array(2,*) ; 5 sec
3: ss = data_array(2,*) ; gsc clean
2: begin ; maccs clean
tf1 = data_array(0,*)
ss = data_array(3,*)
bx1 = long(data_array(04:07,*),0,1,y(type))
bx2 = long(data_array(08:11,*),0,1,y(type))
by1 = long(data_array(12:15,*),0,1,y(type))
by2 = long(data_array(16:19,*),0,1,y(type))
bz1 = long(data_array(20:23,*),0,1,y(type))
bz2 = long(data_array(24:27,*),0,1,y(type))
end
1: begin ; gsc raw
dd = fix (data_array( 0: 1,*),0,1,y(type))
ss = data_array(4,*)
end
0: begin ; maccs raw
header= data_array( 0,*)
dd = fix (data_array( 1: 2,*),0,1,y(type))
ss = data_array(5,*)
bx2 = long(data_array(18:21,*),0,1,y(type))
print, bx2(10)
by2 = long(data_array(22:25,*),0,1,y(type))
print, by2(10)
bz2 = long(data_array(26:29,*),0,1,y(type))
print, bz2(10)
tf1 = data_array(30,*)
tf2 = data_array(31,*)
block = data_array(32,*)
status= data_array(33,*)
house = data_array(34,*)
end
endcase
time = hh*3600. + mm*60. + ss
if keyword_set(int) then begin
good = where(bx1 ne 32767000,count)
if count gt 0 then bx1=interpol(bx1(good),time(good),time)
good = where(by1 ne 32767000,count)
if count gt 0 then bx1=interpol(by1(good),time(good),time)
good = where(bz1 ne 32767000,count)
if count gt 0 then bx1=interpol(bz1(good),time(good),time)
endif
if keyword_set(nt) then begin
bx1 = bx1/1000.
by1 = by1/1000.
bz1 = bz1/1000.
if type eq 2 then begin
bx2 = bx2/1000.
by2 = by2/1000.
bz2 = bz2/1000.
endif
endif
if keyword_set(diag) then print,bx1(0:10),format='(11f7.1)'
end
