#  ADDWCS --- Task to add WCS information to solar images

include <imhdr.h>
define	MAXKEY	80		# Maximum length of value

procedure t_addwcs ()

char	image[SZ_FNAME]

pointer	im

pointer	immap()

begin

#     This program adds World Coordinate System (WCS) information to a
#     solar image.  It attempts to identify, on the basis of keywords in
#     the header, where the image comes from, and then to interpret the
#     image radius and center.
#
#     Author:
#
#     Stephen Walton
#     Department of Physics and Astronomy
#     California State University Northridge 
#     18111 Nordhoff St.
#     Northridge, CA 91330-8268
#     stephen.walton@csun.edu
#
#     Version 1.0.
#

	call clgstr ("image", image, SZ_FNAME)
	im = immap (image, READ_WRITE, 0)

	call addwcs (im)
	call imunmap (im)
end

procedure addwcs(im)
pointer	im
int	imaccf(), imgstr()
bool	streq()
char	result[MAXKEY]
int	junk
begin

	
#	This subroutine adds World Coordinate System (WCS) information to
#	the IRAF file open on unit IM.  The image type is identified
#	first on the basis of the keywords in the header, and then the WCS
#	is added.

#	First, attempt to identify the observatory.  Because of SFO
#	similarities to GONG, we have to handle it specially first.

	if (imaccf(im, "FNDLMBXC") == YES && imaccf(im, "INSTRUME") == YES) {
	    junk = imgstr(im, "INSTRUME", result, MAXKEY)
	    if (streq(result, "CFDT") || streq(result, "CFDT2")) {
	        call printf("Adding WCS to SFO image\n")
	        call sfo_addwcs(im)
	        return
	    }
	}
	if (imaccf(im, "FNDLMBXC") == YES) {
	    call printf("Adding WCS to GONG image\n")
	    call gong_addwcs(im)
	} else if (imaccf(im, "E_XCEN") == YES) {
	    call printf("Adding WCS to KPVT image\n")
	    call kpvt_addwcs(im)
	} else if (imaccf(im, "CENTER_X") == YES) {
	    call printf("Adding WCS to SOHO image\n")
	    call soho_addwcs(im)
	} else {
	    call eprintf("Cannot identify type of image\n")
	}

end

# SFO_ADDWCS --- Add WCS geometry information to an SFO image.

procedure sfo_addwcs(im)
pointer	im
real	xc			# X coordinate of image center
real	yc			# Y coordinate of image center
real	xr			# X radius in pixels
real	yr			# Y radius in pixels
real	fixcfdt			# Tangent of skew angle of x axis
real	radius			# Solar radius in arc seconds
real	pangle			# Solar position angle
real	theta			# Skew angle of x axis
real	ref[2]			# Reference coordinates
real	w[2]			# World coordinates of references coords
real	cd[2, 2]		# CD matrix
real	ltm[2, 2]		# Lterm matrix
real	ltv[2]			# Lterm vector
pointer	mw			# MWCS handle
pointer	mw_openim()
real	imgetr()
int	imaccf()

errchk	mw_newsystem
errchk	imgetr

begin

	iferr {
	    xc = imgetr(im, "FNDLMBXC")
	    yc = imgetr(im, "FNDLMBYC")
	    xr = imgetr(im, "FNDLMBMI")
	    yr = imgetr(im, "FNDLMBMA")
	} then {
	    call eprintf("Did not find one or more required keywords.\n")
	    return
	}
	iferr (fixcfdt = imgetr(im, "FIXCFDT")) fixcfdt = 0.0
	ref[1] = xc
	ref[2] = yc
	w[1] = 0.0
	w[2] = 0.0

#	Check for the value of the solar P angle and add it if not present.

	if (imaccf(im, "SOLAR_P0") == NO) {
	    call addephem (im, radius, pangle)
	    pangle = pangle*acos(-1.0)/180.0
	} else {
	    pangle = imgetr (im, "SOLAR_P0")
	    pangle = pangle*acos(-1.0)/180.0
	}
	theta = atan2(fixcfdt,1.0)

#	The transformation matrix below assumes we are dealing with an
#	unmodified CFDT image.  In such an image, the lower left corner
#	(first datum in image) is the NW corner.  The image is aligned
#	horizontally on geographic North-South and East-West, so it
#	is rotated by the solar P angle.

	cd[1, 1] = -sin(theta)*cos(pangle)/yr - cos(theta)*sin(pangle)/xr
	cd[2, 1] = -cos(pangle)/yr
	cd[1, 2] = sin(theta)*sin(pangle)/yr - cos(theta)*cos(pangle)/xr
	cd[2, 2] = sin(pangle)/yr

	mw = mw_openim (im)
	iferr (call mw_newsystem (mw, "world", 2))
	    ;
	else {
	    call mw_swtype (mw, 1, 1, "linear", "")
	    call mw_swtype (mw, 2, 1, "linear", "")
	}
	ltv[1] = 0.
	ltv[2] = 0.
	ltm[1,1] = 1.
	ltm[1,2] = 0.
	ltm[2,1] = 0.
	ltm[2,2] = 1.
	call mw_sltermr (mw, ltm, ltv, 2)
	call mw_swtermr (mw, ref, w, cd, 2)
	call mw_saveim (mw, im)
	call mw_close (mw)
end

# GONG_ADDWCS - Add WCS information to a GONG image

procedure gong_addwcs(im)
pointer im
real	xc			# X coordinate of image center
real	yc			# Y coordinate of image center
real	xr			# X radius in pixels
real	yr			# Y radius in pixels
real	angle			# Rotation angle of ellipse
real	pi
real	ref[2]			# Reference coordinates
real	w[2]			# World coordinates of references coords
real	cd[2, 2]		# CD matrix
real	ltm[2, 2]		# Lterm matrix
real	ltv[2]			# Lterm vector
pointer	mw			# MWCS handle
pointer	mw_openim()
real	imgetr()

errchk	mw_newsystem
errchk	imgetr
begin
	iferr {
	    xc = imgetr(im, "FNDLMBXC")
	    yc = imgetr(im, "FNDLMBYC")
	    xr = imgetr(im, "FNDLMBMI")
	    yr = imgetr(im, "FNDLMBMA")
	} then {
	    call eprintf("Did not find one or more required keywords.\n")
	    return
	}
	iferr (angle = imgetr(im, "FNDLMBAN")) angle = 0.0
	pi = 4.0*atan(1.0)
	angle = angle*pi/180.0
	ref[1] = xc
	ref[2] = yc
	w[1] = 0.0
	w[2] = 0.0
	
	# The transformation matrix below is actually the product of three
	# matrices, the first of which rotates the system by an angle FNDLMBAN
	# to produce a system with the X axis in the FNDLMBMI direction
	# and the Y axis in the FNDLMBMA direction.  Then, the image is scaled
	# by 1/FNDLMBMI along X and 1/FNDLMBMA along Y, and the image rotated
	# back by -angle
	cd[1,1] = cos(angle)**2/xr + sin(angle)**2/yr
	cd[1,2] = cos(angle)*sin(angle)*(1.0/xr - 1.0/yr)
	cd[2,1] = cd[1,2]
	cd[2,2] = cos(angle)**2/yr + sin(angle)**2/xr

	mw = mw_openim (im)
	iferr (call mw_newsystem (mw, "world", 2))
	    ;
	else {
	    call mw_swtype (mw, 1, 1, "linear", "")
	    call mw_swtype (mw, 2, 1, "linear", "")
	}
	ltv[1] = 0.
	ltv[2] = 0.
	ltm[1,1] = 1.
	ltm[1,2] = 0.
	ltm[2,1] = 0.
	ltm[2,2] = 1.
	call mw_sltermr (mw, ltm, ltv, 2)
	call mw_swtermr (mw, ref, w, cd, 2)
	call mw_saveim (mw, im)
	call mw_close (mw)
end

procedure kpvt_addwcs(im)
pointer im
real	xc			# X coordinate of image center
real	yc			# Y coordinate of image center
real	xr			# X radius in pixels
real	yr			# Y radius in pixels
real	ref[2]			# Reference coordinates
real	w[2]			# World coordinates of references coords
real	cd[2, 2]		# CD matrix
real	ltm[2, 2]		# Lterm matrix
real	ltv[2]			# Lterm vector
pointer	mw			# MWCS handle
pointer	mw_openim()
real	imgetr()

begin
	iferr {
	    xc = imgetr(im, "E_XCEN")
	    yc = imgetr(im, "E_YCEN")
	    xr = imgetr(im, "E_XSMD")
	    yr = imgetr(im, "E_YSMD")
	} then {
	    call eprintf("Did not find one or more required keywords.\n")
	    return
	}
	ref[1] = xc
	ref[2] = yc
	w[1] = 0.0
	w[2] = 0.0
	
	# A very simple transformation.  Scale by 1/r on both axes.
	
	cd[1,1] = 1/xr
	cd[1,2] = 0
	cd[2,1] = 0
	cd[2,2] = 1/yr

	mw = mw_openim (im)
	iferr (call mw_newsystem (mw, "world", 2))
	    ;
	else {
	    call mw_swtype (mw, 1, 1, "linear", "")
	    call mw_swtype (mw, 2, 1, "linear", "")
	}
	ltv[1] = 0.
	ltv[2] = 0.
	ltm[1,1] = 1.
	ltm[1,2] = 0.
	ltm[2,1] = 0.
	ltm[2,2] = 1.
	call mw_sltermr (mw, ltm, ltv, 2)
	call mw_swtermr (mw, ref, w, cd, 2)
	call mw_saveim (mw, im)
	call mw_close (mw)
end

procedure soho_addwcs(im)
pointer im
real	xc			# X coordinate of image center
real	yc			# Y coordinate of image center
real	xr			# X radius in pixels
real	yr			# Y radius in pixels
real	ref[2]			# Reference coordinates
real	w[2]			# World coordinates of references coords
real	cd[2, 2]		# CD matrix
real	ltm[2, 2]		# Lterm matrix
real	ltv[2]			# Lterm vector
pointer	mw			# MWCS handle
pointer	mw_openim()
real	imgetr()

begin
	iferr {
	    xc = imgetr(im, "CENTER_X")
	    yc = imgetr(im, "CENTER_Y")
	    xr = imgetr(im, "R_SUN")
	    yr = xr
	} then {
	    call eprintf("Did not find one or more required keywords.\n")
	    return
	}
	ref[1] = xc
	ref[2] = yc
	w[1] = 0.0
	w[2] = 0.0
	
	# A very simple transformation.  Scale by 1/r on both axes.
	
	cd[1,1] = 1/xr
	cd[1,2] = 0
	cd[2,1] = 0
	cd[2,2] = 1/yr

	mw = mw_openim (im)
	iferr (call mw_newsystem (mw, "world", 2))
	    ;
	else {
	    call mw_swtype (mw, 1, 1, "linear", "")
	    call mw_swtype (mw, 2, 1, "linear", "")
	}
	ltv[1] = 0.
	ltv[2] = 0.
	ltm[1,1] = 1.
	ltm[1,2] = 0.
	ltm[2,1] = 0.
	ltm[2,2] = 1.
	call mw_sltermr (mw, ltm, ltv, 2)
	call mw_swtermr (mw, ref, w, cd, 2)
	call mw_saveim (mw, im)
	call mw_close (mw)
end