C+
C NAME:
C	smei_frm_fts_axis
C PURPOSE:
C	Process header entries for optical axis
C CALLING SEQUENCE:
	subroutine smei_frm_fts_axis(iU, hdr)
C INPUTS:
C	iU		integer			logical unit number of open Fits file
C	hdr(*)		double precision	SMEI header array
C OUTPUTS:
C INCLUDE:
	include		'smei_frm_hdr.h'
C CALLS:
C	smei_hdr_quaternion, quaternion_rotate_xyz
C	FTUKYD, say_fts
C PROCEDURE:
C	Header entries for RA, dec and PA of optical axis are written
C	into the open Fits file.
C MODIFICATION HISTORY:
C	MAR-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		iU
	    double precision	hdr(*)

	character	cSay*17		/'smei_frm_fts_axis'/
	character	cSeverity	/'E'/

	double precision qq(4)
	double precision rx
	double precision ry
	double precision rz
	double precision px
	double precision py
	double precision pz
	double precision phi
	double precision lat
	double precision pa

	double precision datan2d
	double precision dasind

	call smei_hdr_quaternion(hdr,2,qq)	! qq is cam-to-sky quaternion

	rx = 0.0d0				! Unit vector along optical axis in cam frame
	ry = 0.0d0
	rz = 1.0d0

	px = 1.0d0				! Unit vector along long dimension of fov
	py = 0.0d0
	pz = 0.0d0

	call quaternion_rotate_xyz(qq,px,py,pz)
	call quaternion_rotate_xyz(qq,rx,ry,rz) ! Rotate to equatorial sky coordinates

	! RA, Dec of optical axis (z-axis), PA of long dim of fov (x-axis) 
	! Replaces SMEI__HDR_RA, SMEI__HDR_DEC and SMEI__HDR_ANGLE

	! The position angles is an angle measured counter-clockwise relative to ecliptic
	! north as viewed from outside the unit sphere. For

	phi = dmod(datan2d(ry,rx)+360d0,360d0)! RA optical axis in range [0,360]
	lat = dasind(max(-1d0,min(rz,1d0)))	! Dec optical axis
	pa  = datan2d( ry*px-rx*py,pz)	! Position angle with equatorial north

	istat = 0

	call FTUKYD(iU,'RA'   , phi, 7, 'ra of optical axis', istat)
	if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	call FTUKYD(iU,'DEC'  , lat, 7, 'dec of optical axis', istat)
	if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	call FTUKYD(iU,'ANGLE', pa, 7, 'position angle fov', istat)
	if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	return
	end