PRO DOP_FITS,sph,outpath,measurement,line_name,rok,mesic,den,tfirst,tlast,drot,P0,ydiam,speed,progver
;
;Stores a Dopplergram FITS file and a text file with the header.
;INPUTS: sph - Dopplergram array to be stored
;        outpath - output directory (string)
;        measurement - measurement name incl. date (string)
;        line_name - code of spectral line or continuum (string)
;        tfirst,tlast - UT in seconds (float)
;        drot - rotation angle (degrees)
;	 P0 - position angle of solar axis (degrees)
;	 ydiam - solar diameter in pixels
;        speed - scanning speed (string)
;        progver - program and version name
;
; Copyright M. Klvana and M. Sobotka, ASU AVCR Ondrejov, Czech Rep., 2009
; Updated by J. Aboudarham, Obs. Paris, France, 2010:
;	* Added keywords PHYSPARA, SOLAR_P0,CONTACT, WAVEUNIT, CENTER_X and CENTER_Y
;	* Changed the names of the following keywords:
;		- DATE-OBS2 to DATE-END
;		- WAVELENG to WAVELNTH
;	in order to comply with FITS recommendations
; (see: http://bass2000.bagn.obs-mip.fr/Tarbes/spip.php?rubrique22
;	http://heasarc.gsfc.nasa.gov/docs/fcg/standard_dict.html
;	http://tdc-www.harvard.edu/software/wcstools/fitsy2k.html )
;

stcode= ['HA','CN','K3','K1','DOP']
stwave= [6562.8,6558.7,3933.7,3932.3,6562.8]
w=where(stcode eq line_name)
wave=stwave(w(0))

if mesic LT 10 then smes='0'+strtrim(mesic,2) else smes=strtrim(mesic,2)
if den LT 10 then sden='0'+strtrim(den,2) else sden=strtrim(den,2)
date=strtrim(rok,2)+'-'+smes+'-'+sden

ss=round(tfirst mod 60)
mm=fix(tfirst/60.)
hh=mm/60
mm=mm mod 60
if hh LT 10 then shh='0'+strtrim(hh,2) else shh=strtrim(hh,2)
if mm LT 10 then smm='0'+strtrim(mm,2) else smm=strtrim(mm,2)
if ss LT 10 then sss='0'+strtrim(ss,2) else sss=strtrim(ss,2)
time1='T'+shh+':'+smm+':'+sss
ss=round(tlast mod 60)
mm=fix(tlast/60.)
hh=mm/60
mm=mm mod 60
if hh LT 10 then shh='0'+strtrim(hh,2) else shh=strtrim(hh,2)
if mm LT 10 then smm='0'+strtrim(mm,2) else smm=strtrim(mm,2)
if ss LT 10 then sss='0'+strtrim(ss,2) else sss=strtrim(ss,2)
time2='T'+shh+':'+smm+':'+sss

mkhdr,head,sph          ;make the basic FITS header

n1=sxpar(head,'NAXIS1')	;recover the dimensions of the image to define the center
n2=sxpar(head,'NAXIS2')

sxaddpar,head,'DATE-OBS',date+time1,'Observation date, UT start of obs',AFTER='DATE'
sxaddpar,head,'DATE-END',date+time2,'Observation date, UT end of obs',AFTER='DATE-OBS'
sxaddpar,head,'ORIGIN','Coimbra University','Responsible institution',AFTER='DATE-END'
sxaddpar,head,'TELESCOP','Spectroheliograph','Name of telescope',AFTER='ORIGIN'
sxaddpar,head,'CONTACT','A. Garcia','Contact person at the institution',AFTER='TELESCOP'
sxaddpar,head,'BUNIT','m/s','Physical unit',AFTER='CONTACT'
sxaddpar,head,'PHYSPARA','VEL_LONG','Physical parameter represented in the data',AFTER='BUNIT'
sxaddpar,head,'CTYPE1','x-coordinate pix','Solar East-West direction - axis X',AFTER='PHYSPARA'
sxaddpar,head,'CTYPE2','y-coordinate pix','Solar South-North direction - axis Y',AFTER='CTYPE1'
sxaddpar,head,'CDELT1',2.2,'Scale in X arcsec per pixel',AFTER='CTYPE2'
sxaddpar,head,'CDELT2',2.2,'Scale in Y arcsec per pixel',AFTER='CDELT1'
sxaddpar,head,'CROT1',drot,'Angle (deg) used to rotate to East-West',AFTER='CDELT2'
sxaddpar,head,'SOLAR_P0',P0,'Position angle (deg) of solar axis (P0 angle)',AFTER='CROT1'
sxaddpar,head,'WAVELNTH',wave,'Dopplergram wavelength in Angstroems',AFTER='SOLAR_P0'
sxaddpar,head,'WAVEUNIT',-10,'Unit of wavelength in power of ten of meter',AFTER='WAVELNTH'
sxaddpar,head,'SPEED',speed,'Scanning speed',AFTER='WAVEUNIT'
sxaddpar,head,'CENTER_X',n1/2,'Position of the center along X-axis in pixels',AFTER='SPEED'
sxaddpar,head,'CENTER_Y',n2/2,'Position of the center along Y-axis in pixels',AFTER='CENTER_X'
sxaddpar,head,'SOLAR_R',ydiam/2,'Solar disc radius in pixels',AFTER='CENTER_Y'
sxaddpar,head,'COMMENT','   Program '+progver,AFTER='SOLAR_R'

writefits,outpath+measurement+line_name+'.fits',sph,head
print,' SAVED  ',outpath+measurement+line_name+'.fits'
get_lun,unt
openw,unt,outpath+measurement+line_name+'.txt'
printf,unt,head
free_lun,unt

END