pro rdhmifits,directory

;  procedure input the inversion results from HMI, save
;  in a common area, including UT

;  INPUTS:	directory	= AZAM data directory

;  common area for outputs
common com_invertresults, $
    fld, psi, azm, eta0, dmp, dop, bzero, b1mu, cen1, alpha, $
    cct, chisq, occt, fld_er, psi_er, azm_er, $
    cen1_er, alpha_er, xx, yy, slat, slong, $
    yr, mnth, dy, hr, minit, sec, cdelt1, cdelt2,rconv


;  bypass reading data etc, already saved
;goto,jump44

;  HMI inversion data
;  establish relevant part of file name for each parameter
strhmiparams = [ $
	'alpha_mag.fits', 'azimuth.fits', 'chisq.fits', 'damping.fits', $
	'dop_width.fits', 'eta_0.fits', 'field.fits', 'inclination.fits', $
	'src_continuum.fits', 'src_grad.fits', 'vlos_mag.fits', $
	'.alpha_err.fits','.azimuth_err.fits','.field_err.fits', $
	'.inclination_err.fits','.vlos_err.fits' $
	]
strvarname = [ $
	'alpha', 'azm', 'chisq', 'dmp', $
	'dop', 'eta0', 'fld', 'psi', $
	'bzero', 'b1mu', 'cen1', $
	'alpha_er','azm_er','fld_er', $
	'psi_er','cen1_er' $
	]

;  number of inversion parameters
ninvpar = n_elements(strhmiparams)

; Append directory path with '/'.
if(strmid(directory, strlen(directory)-1, 1) ne '/') then begin
	directory = directory + '/'
endif


;  loop over parameters, reading in the FITS files
for kk = 0,ninvpar-1 do begin
	derf = file_search(directory + '*' + strhmiparams(kk))
	filnam = derf
	tdat = readfits(filnam,hdr)
;  eliminate NaNs, changing them to zero values
	whr = where(tdat ne tdat,ncount)
	if (ncount ne 0) then tdat(whr) = 0.
;  assign names to the variable
	void = execute(strvarname[kk] + '= tdat')
;  for first variable, assign RCONV, unity for fitted points
	if kk eq 0 then rconv = alpha
endfor

;;  !!!!temporarily use the computed continuum for continuum
;cct = bzero+b1mu
;occt = cct

;;  read in the extreme raw intensity images, average them
;;  search data directory for I0 and I5
;goto, jump55
;derfi0 = file_search(directory + '*' + 'I0' + '*' + '.fits')
;;  read in I0 continuum intensity
;cont0 = fitsio_read_image(derfi0(0),hdrc0)
;whr = where(cont0 ne cont0)
;cont0(whr) = 0.
;;  read in I5 continuum intensity
;derfi5 = file_search(directory + '*' + 'I5' + '*' + '.fits')
;cont5 = fitsio_read_image(derfi5(0),hdrc5)
;whr = where(cont5 ne cont5)
;cont5(whr) = 0.
;;  average the images
;occt = float(0.5*(cont0 + cont5))
;jump55:

;  read in the continuum intensity from observed data
derfi0 = file_search(directory + 'continuum.fits')
cont0 = fitsio_read_image(derfi0(0),hdrc0)
whr = where(cont0 ne cont0)
cont0(whr) = 0.
occt = float(cont0)
cct = occt

;  constrain field > 0
fld = fld>0.

;  rotate azimuth angle so that zero azimuth is solar west.
;  The inversion code places zero azimuth along pixel column,
;  increasing counter-clockwise
azm = azm + 90.



;  save HMI FITS header in the  AZAM directory as an
;  IDL save-file
header = hdr
save,filename=directory+'HMI_FITSheader.save',header

;  get the time information from the FITS file header
timez = sxpar(hdr,'T_OBS')
yr=fix(strmid(timez,0,4))
mnth=fix(strmid(timez,5,2))
dy=fix(strmid(timez,8,2))
hr=fix(strmid(timez,11,2))
minit=fix(strmid(timez,14,2))
sec=float(strmid(timez,17,6))

;  find the sun center point from the header values
;  note:  first pixel is 1, not 0! Remove this offset for IDL
xctr = sxpar(hdr,'CRPIX1') - 1.
yctr = sxpar(hdr,'CRPIX2') - 1.

;  find the rotation parameter, presume CLOCKWISE from solar N, degrees
angrot = sxpar(hdr,'CROTA2')
;  if the rotation parameter is close to 180 degrees, then the HMI
;  images will be rotated by exactly 180 degrees because we will eventually
;  need to specify the azimuth for the original data presentation, not
;  the rotated one.  Must do the rotation otherwise the field angles
;  will not be presented to the user correctly
if(abs(angrot-180.) le 4.0) then begin
	for kk=0,ninvpar-1 do begin
		void = execute(strvarname[kk] + '= rotate(' + strvarname[kk] + ',2)' )
	endfor 
endif else begin
	print,'CROTA2 = ',angrot,' > 4 deg. Cannot display in AZAM correctly'
	stop
endelse
;  rotate the continuum image
occt = rotate(occt,2)
cct = occt

;  with 180 deg rotation, now must find new sun center positions
;  get some dimensions
npix = n_elements(fld)
sz = size(fld)
nxx = sz(1)
nyy = sz(2)
xctr = float(nxx)-1-xctr
yctr = float(nyy)-1-yctr

;  get arrays for x,y positions, arcsec
cdelt1 = sxpar(hdr,'CDELT1')
xx = fld & xx(*,*) = 0.
cdelt2 = sxpar(hdr,'CDELT2')
for ii = 0,nxx-1 do xx(ii,*) = (float(ii) - xctr)*cdelt1
yy = fld & yy(*,*) = 0.
for jj = 0,nyy-1 do yy(*,jj) = (float(jj) - yctr)*cdelt2



;  find the solar latitude,longitude from XX, YY
;  xy2lonlat calculates its own ephemeris values based on times
;  make a loop to try to speed up this really slow calculation
slong = fltarr(nxx,nyy) & slat = fltarr(nxx,nyy)
;  process by column
arcvec = fltarr(2,nxx)
time1d = strarr(nxx)
time1d(*) = timez
print,' NOTICE: beginning computation of solar latitude,longitude.'
print,' The SolarSoft routine xy2lonlat is VERY SLOW.....have patience!'
for jj = 0,nyy-1 do begin
;  calculate modulus, and print if 100
	jjmod = jj mod 100
	if jjmod eq 0 then print,'starting lat,long computation for row ',jj
	arcvec(0,*) = xx(*,jj)
	arcvec(1,*) = yy(*,jj)
	lonlat = xy2lonlat(arcvec,time1d)
	slong(*,jj) = lonlat(0,*)
	slat(*,jj) = lonlat(1,*)
endfor

;save,filename='rdhmifits.save', $
;	fld, psi, azm, eta0, dmp, dop, bzero, b1mu, cen1, alpha, $
;    cct, chisq, occt, fld_er, psi_er, azm_er, $
;    cen1_er, alpha_er, xx, yy, slat, slong, $
;    yr, mnth, dy, hr, minit, sec, cdelt1, cdelt2,rconv

;jump44: restore,'rdhmifits.save'

return
end