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

;  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[0]
;	tdat = readfits(filnam,hdr)

;  use new IDL HMI read package
read_sdo,derf,hdr,tdat

;  eliminate NaNs, changing them to zero values
	whrnan = where(tdat ne tdat,ncount)
	if (ncount ne 0) then tdat(whrnan) = 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

;  read in the continuum intensity from observed data
derfi0 = file_search(directory + 'continuum.fits')

;cont0 = fitsio_read_image(string(derfi0(0)),hdrc0)
;cont0 = fitsio_read_image(derfi0(0),hdrc0)
;  use new IDL HMI read package
read_sdo,derfi0,hdrc0,cont0

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
;  ???NOTE: read_sdo saves header as a structure, so sxpar does not work
;  below we extract the desired header values from the structure instead)
;  Therefore, header is not saved here
;header = hdr
;save,filename=directory+'HMI_FITSheader.save',header

;  get the time information from the FITS file header
;timez = sxpar(hdr,'T_OBS')
timez = 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.
xctr = hdr.CRPIX1 - 1.
yctr = hdr.CRPIX2 - 1.

;  find the rotation parameter, presume CLOCKWISE from solar N, degrees
;angrot = sxpar(hdr,'CROTA2')
angrot = 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')
cdelt1 = hdr.CDELT1
xx = fld & xx(*,*) = 0.
;cdelt2 = sxpar(hdr,'CDELT2')
cdelt2 = 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
utimed = double(float(hr) + float(minit)/60. + sec/3600.)
;  expand time arrays to full array size, not just per slit position
iyr = intarr(nxx,nyy) & imnth = intarr(nxx,nyy) & idy = intarr(nxx,nyy)
ihr = intarr(nxx,nyy) & iminit = intarr(nxx,nyy)
secc = fltarr(nxx,nyy) & uutimed =fltarr(nxx,nyy)
iyr(*,*) = yr 
imnth(*,*) = mnth
idy(*,*) = dy
ihr(*,*) = hr
iminit(*,*) = minit
secc(*,*) = sec
uutimed = utimed
azam_sunframe, iyr, imnth, idy, uutimed, xx, yy, $
	slong, slat, valuemu, $
    ras, dcs, gsdt, b0ang, peeang, ctrlong, r0
;  convert longitude,latitude from radians to degrees
	slong = float(slong/!dtor)
	slat = float(slat/!dtor)

return
end