;set of IDL routines to read VIRGO data files
;
;	needs IDL FITS library
;	if you encounter any problems please let me know!
;	Version 1.1: prepared 9-May-2003 by Cfrohlich@pmodwrc.ch


;soho mission day
function ymd2soho_md, yymmdd
	y=long(yymmdd)/10000L
	m=(yymmdd-y*10000L)/100L
	d=yymmdd-y*10000L-m*100L
	if y gt 50 then y=y+1900L else y=y+2000L
	return, julday(m,d,y)-julday(12,1,1995)
end
;**************************************************************

pro read_level1P,filename,day,value,flag
;Procedure to read daily PMO6V FITS-Files
;returns 1440 values of TSI and flag
  d_temp=readfits(filename,header,exten=1,/silent)
  read_header,header,'NAXIS2  ',tempd
  if strmid(filename,1,2,/revers) eq '00' then $
  	read_header,header,'TNULL1  ',notval else $
  	read_header_str,header,'TNULL1  ',notval
  n=long(tempd)
  date=long(strmid(filename,11,6,/revers))
  day=dindgen(n)/tempd+double(ymd2soho_md(date))+22.d0/86400.d0
  value=double(ftget(header,d_temp,'PMO6-VA '))
  flag=ftget(header,d_temp,'DQ PMO6-VA')
  inval=where(value le notval,cnt)
  if cnt gt 0 then value(inval)=!values.D_NaN
end
;**************************************************************

pro read_level1D,filename,day,value
;Procedure to read daily DIARAD FITS-Files
;returns 1440 values of TSI
  d_temp=readfits(filename,header,exten=1,/silent)
  read_header,header,'NAXIS2  ',tempd
  read_header_str,header,'TNULL1  ',notval
  n=long(tempd)
  date=long(strmid(filename,11,6,/revers))
  day=double(ymd2soho_md(date))+dindgen(n)/tempd+87.d0/86400.d0
  value=double(ftget(header,d_temp,'SOLAR_CT'))
  inval=where(value le notval,cnt)
  if cnt gt 0 then value(inval)=!values.D_NaN
end
;**************************************************************

pro read_level1S,filename,day,value,flag
;Procedure to read daily SPM FITS-Files
;returns 1440 values of red, green and blue channel and flag
  d_temp=readfits(filename,header,exten=1,/silent)
  read_header,header,'NAXIS2  ',tempd
  read_header_str,header,'TNULL1  ',notval
  n=long(tempd)
  date=long(strmid(filename,11,6,/revers))
  day=double(ymd2soho_md(date))+dindgen(n)/tempd+32.d0/86400.d0
  value=dblarr(3,n)
  value[0,*]=double(ftget(header,d_temp,'RED     '))
  value[1,*]=double(ftget(header,d_temp,'GREEN   '))
  value[2,*]=double(ftget(header,d_temp,'BLUE    '))
  flag=ftget(header,d_temp,'DQ1     ')
  inval=where(value ge notval,cnt)
  if cnt gt 0 then value(inval)=!values.D_NaN
end
;**************************************************************

pro read_level1L,filename,day,value,flag,diameter
;Procedure to read daily LOI FITS-Files
;returns 1440 values of 12 inner pixel and 4 guiding pixel and 2 flags for the two groups
;		and of 2 diameters
  d_temp=readfits(filename,header,exten=1,/silent)
  read_header,header,'NAXIS2  ',tempd
  read_header_str,header,'TNULL1  ',notval
  n=long(tempd)
  date=long(strmid(filename,11,6,/revers))
  day=double(ymd2soho_md(date))+dindgen(n)/tempd+32.d0/86400.d0
  value=dblarr(16,n)
  flag=strarr(2,n)
  for i=0,15 do begin
  		if i lt 9 then label='PIXEL'+string(i+1,format='(i1)')+'  ' else $
  				   label='PIXEL'+string(i+1,format='(i2)')+' '
  		value[i,*]=double(ftget(header,d_temp,label))
  endfor
  flag(0,*)=ftget(header,d_temp,'QFLAGSC ')
  flag(1,*)=ftget(header,d_temp,'QFLAGGU ')
  inval=where(value ge notval,cnt)
  if cnt gt 0 then value(inval)=!values.D_NaN
  diameter=dblarr(2,n)
  diameter(0,*)=ftget(header,d_temp,'DIAM_EW ')
  diameter(1,*)=ftget(header,d_temp,'DIAM_NS ')
end
;**************************************************************

;read fits header according to keyword as string
pro read_header_str, header, keyword, var
ind=where(strmid(header,0,8) eq keyword,cnt)
if cnt ne 1 then begin
	print, 'READ_HEADER: no such keyword: ',keyword
	var=!values.F_NaN
	return
endif
comment_pos=strpos(header(ind(0)),'/',9)
if comment_pos gt 0 then $
	var_string=strmid(header(ind(0)),9,comment_pos-10L) $
else var_string=strmid(header(ind(0)),9,24)
var_string=strtrim(var_string,2)
var=float(strmid(var_string,1,strlen(var_string)-2))
return
end
;**************************************************************

;read fits header according to keyword as value
pro read_header, header, keyword, var
ind=where(strmid(header,0,8) eq keyword,cnt)
if cnt ne 1 then begin
	print, 'READ_HEADER: no such keyword: ',keyword
	var=''
	return
endif
comment_pos=strpos(header(ind(0)),'/',9)
if comment_pos gt 0 then $
	var_string=strmid(header(ind(0)),9,comment_pos-10L) $
else var_string=strmid(header(ind(0)),9,24)
reads,var_string,var
return
end
;**************************************************************

pro test_fits

dir='D:/Data/VIRGO/Examples of Level1 FITS Files/'
;change this to the path where you fits files are

;test files used
file_spm=dir+'030101_1.S02'
file_loi=dir+'030101_1.L00'
file_dia=dir+'030101_1.D01'

;in the header of *.P00 the invalid is given as value
;	and in *.P01 as string (same format as the other VIRGO FITS headers)
file_pmo0=dir+'030101_1.P00'
file_pmo1=dir+'030430_1.P01'

read_level1P,file_pmo0,day,value,flag
print,'PMO6V old header:'
help,value
print,day[0:5],value[0:5],flag[0:5]
stop

read_level1P,file_pmo1,day,value,flag
print,'PMO6V new header:'
help,value
print,day[0:5],value[0:5],flag[0:5]
stop

read_level1D,file_DIA,day,value
print,'DIARAD:'
help,value
print,'first 60 values are NaN due to DIARAD calibration'
print,day[0:5],value[0:5]
print,day[60:65],value[60:65]
stop

read_level1S,file_spm,day,value,flag
print,'SPM:'
help,value
print,day[0:5],reform(value[0,0:5]),reform(value[1,0:5]),reform(value[2,0:5]),flag[0:5]
stop

read_level1L,file_LOI,day,value,flag,diameter
print,'LOI:'
help,value
print,day[0:5],reform(value[0,0:5]),reform(value[12,0:5]),reform(flag[0,0:5]),$
	reform(flag[1,0:5]),reform(diameter[0,0:5]),reform(diameter[1,0:5])

stop

end