pro rdmerlinfits,directory,merdat

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

;  INPUTS:	directory	= AZAM data directory
;				merdat	= path name to the Merlin output fits file

;  common area for outputs
common com_merlinresults, $
	fld, psi, azm, eta0, dmp, dop,$
	bzero, b1mu, vmac, cen1, cen2, alpha, $
	delta, cct, chisq, chisq_i, chisq_q, $
	chisq_u, chisq_v, signed_v, niter, occt, $
	rconv, netpol, netpol_full, $
	fld_er, psi_er, azm_er, eta0_er, dmp_er, dop_er, $
	bzero_er, b1mu_er, vmac_er, cen1_er, cen2_er, $
	alpha_er, delta_er,xx,yy,slat,slong, $
	yr,mnth,dy,hr,minit,sec



;  MERLIN data
strr = azam_merlin_fits_read(directory,merdat)


;  setup arrays
fld = strr.field_strength
psi = strr.field_inclination
azm = strr.field_azimuth
eta0 = strr.line_strength
dmp = strr.damping
dop = strr.doppler_width
bzero = strr.source_function
b1mu = strr.source_function_gradient
vmac = strr.macro_turbulence
cen1 = strr.doppler_shift1
cen2 = strr.doppler_shift2
alpha = strr.stray_light_fill_factor
delta = strr.stray_light_shift
cct = strr.continuum_intensity
occt = strr.original_continuum_intensity
chisq = strr.chisq_total
chisq_i = strr.chisq_i
chisq_q = strr.chisq_q
chisq_u = strr.chisq_u
chisq_v = strr.chisq_v
signed_v = strr.stokesv_magnitude
niter = strr.number_of_iterations
rconv = strr.fitting_attribute
netpol = strr.polarization
;mer_netpol_full = dblarr(nmer,nsl)
fld_er = strr.field_strength_error
psi_er = strr.field_inclination_error
azm_er = strr.field_azimuth_error
eta0_er = strr.line_strength_error
dmp_er = strr.damping_error
dop_er = strr.doppler_width_error
bzero_er = strr.source_function_error
b1mu_er = strr.source_function_gradient_error
vmac_er = strr.macro_turbulence_error
cen1_er = strr.doppler_shift1_error
cen2_er = strr.doppler_shift2_error
alpha_er = strr.stray_light_fill_factor_error
delta_er = strr.stray_light_shift_error
xx = strr.x_coordinate
yy = strr.y_coordinate

;  get some dimensions
npix = n_elements(xx)
sz = size(xx)
nxx = sz(1) & nyy = sz(2)

;  get the time information from the FITS file
timez = strr.times.times
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 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 # replicate(1., nyy)
imnth(*,*) = mnth # replicate(1., nyy)
idy(*,*) = dy # replicate(1., nyy)
ihr(*,*) = hr # replicate(1., nyy)
iminit(*,*) = minit # replicate(1., nyy)
secc(*,*) = sec # replicate(1.,nyy)
uutimed = utimed # replicate(1.,nyy)
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)

;;  old slow wrong method
;;  find the solar latitude,longitude from XX, YY
;mdpt = n_elements(hr)/2
;timmm = timez(mdpt)
;;  get a 2-D string array for time to enable more precise
;;  ephemeris calculations throughout maps
;time2d = strarr(nxx,nyy)
;for jj = 0,nyy-1 do time2d(*,jj) = timez
;;  xy2lonlat calculates its own ephemeris values based on times
;arcvec = fltarr(2,npix)
;arcvec(0,*) = xx
;arcvec(1,*) = yy
;lonlat = xy2lonlat(arcvec,time2d)
;slong = fltarr(nxx,nyy) & slat = fltarr(nxx,nyy)
;slong(*,*) = lonlat(0,*)
;slat(*,*) = lonlat(1,*)

return
end