; Takes a complete MDI magnetogram and remaps with corrected Bz pixel values
; based on helioscopic angle
pro bzcorrect, fname

	data=readfits(fname, header)
        x0=fxpar(header,'X0')
        y0=fxpar(header,'Y0')
	xscale=fxpar(header,'XSCALE')
        yscale=fxpar(header,'YSCALE')
        angle=fxpar(header,'P_ANGLE')
        blnk=fxpar(header,'BLANK')
	date=strmid(fxpar(header,'T_START'),0,10)

	blanks=where(data le blnk)

; Adjust mdi pixel values using angle from solar disc center
	x_values=fltarr(n_elements(data(*,0)),n_elements(data(0,*)))
	for i=0, n_elements(data(*,0))-1 do x_values(i,0:n_elements(data(0,*))-1)=dindgen(1,n_elements(data(0,*)))
        y_values=fltarr(n_elements(data(*,0)),n_elements(data(0,*)))
        for i=0, n_elements(data(0,*))-1 do y_values(0:n_elements(data(*,0))-1,i)=dindgen(1,n_elements(data(0,*)))

	arcsec=fltarr(2,n_elements(x_values))
	arcsec(0,*)=(x_values-X0)*XSCALE & arcsec(1,*)=(y_values-Y0)*YSCALE

	lonlat=conv_a2h(arcsec, date, off_limb=off_limb)
	central_angle=fltarr(n_elements(data(*,0)),n_elements(data(0,*)))
	for i=0L, n_elements(central_angle)-1 do central_angle(i)=acos( cos(lonlat(0,i)*0.0174532925) * cos(lonlat(1,i)*0.0174532925) )

	data=data/(sin(acos(central_angle*0.636619772)))

; Correct data for MDI scaling issues

 	hi=where(abs(data) gt 1200.)
	if hi(0) ne -1 then begin
		hi=(abs(data) gt 1200.)*1.9
		lo=(abs(data) le 1200.)*1.45
		data=data*(hi+lo)
	endif else data=data*1.45

; Overwrite corrected fits file
	spawn, 'rm -f '+fname
	writefits, fname, float(data), header, nanvalue=blnk
	
end