# File rvsao/Util/t_bcvcorr.x
# December 8, 1999
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics

# Copyright(c) 1995-1999 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose.  It is provided "as is"
# without express or implied warranty.

# BCVCORR is an IRAF task which computes heliocentric and solar system
# barycentric corrections for spectra


include	<imhdr.h>
include	<fset.h>
include <imhdr.h>
include <imio.h>

define	LEN_USER_AREA	100000
define	MAX_RANGES	50
define	SZ_HKWORD	9
define	SZ_MONTH	3

procedure t_bcvcorr ()

pointer	specim		# image structure for spectrum
char	spectra[SZ_LINE]	# Object spectra list
char	specfile[SZ_PATHNAME]	# Object spectrum file name
char	specdir[SZ_PATHNAME]	# Directory for object spectra
char	specpath[SZ_PATHNAME]	# Object spectrum path name
char	obsname[SZ_LINE]	# Observatory name
char	keyword[SZ_LINE]	# Various keywords
char	keyobs[SZ_LINE]		# Observatory name keyword
char	keyhjd[SZ_LINE]		# Header keyword for Heliocentric Julian Day
char	keyjd[SZ_LINE]		# Header keyword for Julian Day
pointer	speclist	# List of spectrum files
bool	savebcv,savebcv0	# Save velocity correction in data file header
int	mode		# Mode (READ_ONLY, READ_WRITE)
pointer	obs		# pointer to observatory structure
bool	gottime		# True if observation time obtained
bool	gotstart,gotend	# True if start or stop observation times obtained
bool	gotdate		# True if observation date obtained
bool	gotpos		# True if position of observed object obtained
#bool	debug
bool	verbose, debug, clgetb()
int	mm, dd, yyyy, ldir, ndim, clgeti()
double	ra, dec, eq, exp, dindef
double	ut1, ut2, utm, ut, dlong, dlat, dalt
double	hjd		# Heliocentric Julian Date
double	gjd		# Geocentric Julian Date
double	dpi,dtr, lt
double	dgcv, dbcv, dhcv, fbcv, fhcv, tpi
real	rbcv, rhcv
int	lspec
int	strlen(), imtgetim(), imaccess(), strncmp()
double	juldat(), clgetd()
char	specnums[SZ_LINE]
int	mspec_range[3,MAX_RANGES]
int	nmspec0,mspec
char	str[SZ_LINE]
char	uts[32]
string  month   "JanFebMarAprMayJunJulAugSepOctNovDec"

pointer	imtopen()
pointer	immap()
errchk	immap()
double	obsgetd()
pointer	obsvopen()
errchk	obsvopen
int	decode_ranges(),get_next_number()
int	imaccf()

define	newspec_	10
define	bcvcomp_	50
define	endvc_	90

begin
	dpi = 3.1415926535897932d0
	tpi = 2.d0 * dpi
	dtr = dpi / 180.0d0
	gjd = 0.d0
	hjd = 0.d0
	gotdate = FALSE
	gottime = FALSE
	gotstart = FALSE
	gotend = FALSE
	verbose = clgetb ("verbose")

# Spectra for which to compute solar system center velocity corrections
	call clgstr ("spectra",spectra,SZ_PATHNAME)
	lspec = strlen (spectra)
	if (lspec == 0) {
	    ra = clgetd ("ra")
	    dec = clgetd ("dec")
	    eq = clgetd ("equinox")
	    if (ra != dindef && dec != dindef)
		gotpos = TRUE
	    else
		gotpos = FALSE
	    hjd = clgetd ("hjd")
	    gjd = clgetd ("gjd")
	    yyyy = clgeti ("year")
	    mm = clgeti ("month")
	    dd = clgeti ("day")
	    if (yyyy != INDEFI && mm != INDEFI && dd != INDEFI)
		gotdate = TRUE
	    else if (gjd != 0.d0 && gjd != dindef)
		gotdate = TRUE
	    else
		gottime = FALSE
	    ut = clgetd ("ut")
	    if (ut != dindef)
		gottime = TRUE
	    else if (gjd != 0.d0 && gjd != dindef)
		gotdate = TRUE
	    else
		gottime = FALSE
	    specim = NULL
	    savebcv = FALSE
	    go to bcvcomp_
	    }
	speclist = imtopen (spectra)
	call clgstr ("specdir",specpath,SZ_PATHNAME)

# Multispec spectrum numbers (use only first if multiple files)
	call clgstr ("specnum",specnums,SZ_LINE)
	if (specnums[1] == EOS)
	    call strcpy ("0",specnums,SZ_LINE)
	call flush (STDOUT)
	if (decode_ranges (specnums, mspec_range, MAX_RANGES, nmspec0) == ERR){
	    call sprintf (str, SZ_LINE, "T_XCSAO: Illegal multispec list <%s>")
	    call pargstr (specnums)
	    call error (1, str)
	    }

	call clgstr ("specdir",specdir,SZ_PATHNAME)
	debug = clgetb ("debug")
	savebcv0 = clgetb ("savebcv")

# Get next object spectrum file name from the list
newspec_
	savebcv = savebcv0
	dindef = INDEFD
	gjd = 0.d0
	hjd = 0.d0
	if (imtgetim (speclist, specfile, SZ_PATHNAME) == EOF)
	   go to endvc_

	if (get_next_number (mspec_range, mspec) == EOF)
	    mspec = 0

# Check for readability of object spectrum
	ldir = strlen (specdir)
	if (ldir > 0) {
	    if (specdir[ldir] != '/') {
		specdir[ldir+1] = '/'
		specdir[ldir+2] = EOS
		}
	    call strcpy (specdir,specpath,SZ_PATHNAME)
	    call strcat (specfile,specpath,SZ_PATHNAME)
	    }
	else
	    call strcpy (specfile,specpath,SZ_PATHNAME)

	if (imaccess (specpath, READ_ONLY) == NO) {
	    call eprintf ("BCVCORR: cannot read spectrum file %s \n")
		call pargstr (specpath)
	    go to newspec_
	    }

# Load spectrum header
	mode = READ_WRITE
        iferr (specim = immap (specpath, mode, LEN_USER_AREA)) {
            mode = READ_ONLY
	    iferr (specim = immap (specpath, mode, LEN_USER_AREA)) {
		if (specim != NULL) call imunmap (specim)
		call eprintf ("BCVCORR: Cannot read image %s\n")
		call pargstr (specpath)
		specim = ERR
		go to endvc_
		}
	    call eprintf ("BCVCORR: Cannot write to image %s\n")
	    call pargstr (specpath)
	    savebcv = FALSE
	    }
        ndim = IM_NDIM(specim)
	if (debug) {
	    call printf ("%s: %d x %d x %d %d-D image\n")
		call pargstr (specpath)
		call pargi (IM_LEN(specim,1))
		call pargi (IM_LEN(specim,2))
		call pargi (IM_LEN(specim,3))
		call pargi (ndim)
	    }


# Direction for velocity correction
	call clgstr ("keyra",keyword,SZ_LINE)
	ra = dindef
 	call imgdpar (specim, keyword, ra)
	call clgstr ("keydec",keyword,SZ_LINE)
	dec = dindef
	call imgdpar (specim,keyword, dec)
	call clgstr ("keyeqnx",keyword,SZ_LINE)
	eq = 1950.0d0
	call imgdpar (specim,keyword,eq)
	if (ra != dindef && dec != dindef)
	    gotpos = TRUE
	else
	    gotpos = FALSE

# Midtime of observation
	call clgstr ("keyhjd",keyhjd,SZ_LINE)
	call imgdpar (specim, keyhjd, hjd)
	call clgstr ("keyjd",keyjd,SZ_LINE)
	call imgdpar (specim, keyjd, gjd)
	if (gjd == 0.0 || gjd == dindef) {
	    call clgstr ("keymid",keyword,SZ_LINE)
	    ut = dindef
	    call imgdpar (specim,keyword, ut)
	    if (ut != dindef)
		gottime = TRUE
	    else {
		call clgstr ("keyend",keyword,SZ_LINE)
		ut2 = dindef
		call imgdpar (specim, keyword, ut2)
		if (ut2 != dindef)
		    gotend = TRUE
		ut1 = dindef
		call clgstr ("keystart",keyword,SZ_LINE)
		call imgdpar (specim,keyword, ut1)
		if (ut1 != dindef)
		    gotstart = TRUE
		exp = 0.
		call clgstr ("keyexp",keyword,SZ_LINE)
		call imgdpar (specim,keyword, exp)
		if (ut1 != dindef && ut2 != dindef) {
		    ut = (ut1 + ut2) * 0.5
		    gottime = TRUE
		    }
		else if (ut2 != dindef && exp > 0.d0) {
		    ut = ut2 - (exp * 0.5 / 3600.)
		    gottime = TRUE
		    }
		else if (ut1 != dindef && exp > 0.d0) {
		    ut = ut1 + (exp * 0.5 / 3600.)
		    gottime = TRUE
		    }
		else {
		    ut = dindef
		    gottime = FALSE
		    }
		if (debug ) {
		    call printf("UT start: %h, mid: %h, end: %h, exp: %d\n")
			call pargd (ut1)
			call pargd (ut)
			call pargd (ut2)
			call pargd (exp)
		    }
		}

# Date of observation, if midtime is not Julian Date
	    mm = INDEFI
	    dd = INDEFI
	    yyyy = INDEFI
	    call clgstr ("keydate",keyword,SZ_LINE)
	    call imgdate (specim, keyword, mm, dd, yyyy)
	    if (ut == dindef) {
		if (gotend)
		    ut = ut2
		else if (gotstart)
		    ut = ut1
		}

# Print date, time, and pointing direction
	    if (debug) {
		call printf("%04d-%3.3s-%02d %h UT")
		    call pargi (yyyy)
		    call pargstr (month[(mm-1) * SZ_MONTH + 1])
		    call pargi (dd)
		    call pargd (ut)
		if (!gottime && gotend)
		    call printf (" end of obs.\n")
		else if (!gottime && gotstart)
		    call printf (" start of obs.\n")
		else
		    call printf ("\n")
		}
	    if (mm != dindef && dd != dindef && yyyy != INDEFI)
		gotdate = TRUE
	    else
		gotdate = FALSE
	    }
	else {
	    if (debug) {
		call printf("Julian date: %.4f\n")
		    call pargd (gjd)
		}
	    gotdate = TRUE
	    }
	if (verbose) {
	    call printf("RA: %h, Dec: %h  %6.1f\n")
		    call pargd (ra)
		    call pargd (dec)
		    call pargd (eq)
	    }
       	gjd = juldat (mm, dd, yyyy) - 0.5d0
       	gjd = gjd + ut / 24.0d0
#	call juldate (specim, ut, gjd, hjd, debug)

# Location of observatory
bcvcomp_
	obsname[1] = EOS
	dlong = dindef
	dlat = dindef
	dalt = 0.d0
	call clgstr ("obsname",obsname,SZ_LINE)
	call clgstr ("keyobs",keyobs,SZ_LINE)

# Read position from image header
	if (specim!=NULL && specim!=ERR && strncmp (obsname,"file",4)==0) {
 	    call imgspar (specim,keyword,obsname,SZ_LINE)
	    if (obsname[1] != EOS) {
		call clgstr ("keylong",keyword,SZ_LINE)
 		call imgdpar (specim,keyword,dlong)
		call clgstr ("keylat",keyword,SZ_LINE)
 		call imgdpar (specim,keyword,dlat)
		call clgstr ("keyalt",keyword,SZ_LINE)
 		call imgdpar (specim,keyword,dalt)
		}
	    }

# Read position from BCVCORR parameters
	if (dlong == dindef || dlat == dindef) {
	    dlong = clgetd ("obslong")
	    dlat = clgetd ("obslat")
	    dalt = clgetd ("obsalt")
	    }

# Read position from IRAF observatory database using name from file
	if ((dlong == dindef || dlat == dindef) && keyobs[1] != EOS) {
 	    call imgspar (specim,keyobs,obsname,SZ_LINE)
	    if (obsname[1] != EOS) {
        	obs = obsvopen (obsname, NO)
		if (obs != NULL) {
#		    call obslog (obs,"BCVCORR","latitude longitude altitude",					 STDOUT)
		    dlat = obsgetd (obs, "latitude")
		    dlong = obsgetd (obs, "longitude")
		    dalt = obsgetd (obs, "altitude")
		    call obsclose (obs)
		    }
		}
	    }

# Read position from IRAF observatory database using name from parameter list
	if ((dlong == dindef || dlat == dindef) && obsname[1] != EOS) {
            obs = obsvopen (obsname, NO)
	    if (obs != NULL) {
#		call obslog (obs, "BCVCORR", "latitude longitude altitude", STDOUT)
		dlat = obsgetd (obs, "latitude")
		dlong = obsgetd (obs, "longitude")
		dalt = obsgetd (obs, "altitude")
		call obsclose (obs)
		}
	    }

	if (verbose) {
	    call printf("%s lat %h , long %h, alt %.1f\n")
		call pargstr (obsname)
		call pargd (dlat)
		call pargd (dlong)
		call pargd (dalt)
	    }
	call flush (STDOUT)

# Reformat so juldat and heliocvel can use the info
	rhcv = 0.
	rbcv = 0.
	if (gotdate && gottime) {

	# Calculate Julian Date
	    if (yyyy == INDEFI || mm == INDEFI || dd == INDEFI || ut == dindef)
		call jd2ut (gjd, yyyy, mm, dd, ut)
	    else if (gjd == 0.0 || gjd == dindef) {
        	gjd = juldat (mm,dd,yyyy) - 0.5d0
        	gjd = gjd + ut / 24.0d0
		}
	    if (verbose) {
		call printf("Julian date is %.5f at %4d-%3.3s-%02d %h UT")
		    call pargd (gjd)
		    call pargi (yyyy)
		    call pargstr (month[(mm-1) * SZ_MONTH + 1])
		    call pargi (dd)
		    call pargd (ut)
		if (!gottime && gotend)
		    call printf (" at end of obs.\n")
		else if (!gottime && gotstart)
		    call printf (" at start of obs.\n")
		else
		    call printf ("\n")
		call flush (STDOUT)
		}

	# Compute the Julian Date at the time when the light reached the Sun
	    if ((hjd == 0.0 || hjd == dindef) &&
		ra != dindef && dec != dindef) {
		call jd2hjd (ra, dec, gjd, lt, hjd)
		}
	    if (verbose ) {
		call printf("Object at ra %.3h dec %.2h eq %.1d\n")
		    call pargd (ra)
		    call pargd (dec)
		    call pargd (eq)
		call printf("Heliocentric Julian date: %.5f\n")
		    call pargd (hjd)
		call flush (STDOUT)
		}

	# Compute radial velocity corrections
	    call bcv (gjd, dlong,dlat,dalt, ra,dec,eq, dbcv, dhcv, dgcv)

	    if (verbose) {
		call printf ("gbcvel = %.4f  ghcvel = %.4f  geovel = %.4f\n")
		    call pargd (dbcv)
		    call pargd (dhcv)
		    call pargd (dgcv)
		call flush (STDOUT)
		}

	# Add up both corrections
	    rbcv = dbcv + dgcv
	    rhcv = dhcv + dgcv
	    if (verbose) {
		call printf ("bcv = %.4f  hcv = %.4f computed\n")
		    call pargr (rbcv)
		    call pargr (rhcv)
		fbcv = dindef
		fhcv = dindef
		if (specim != NULL && specim != ERR) {
		    call imgdpar (specim, "BCV", fbcv)
		    call imgdpar (specim, "HCV", fhcv)
		    call printf ("bcv = %.4f  hcv = %.4f from file\n")
			call pargd (fbcv)
			call pargd (fhcv)
		    }
		call flush (STDOUT)
		}

# Check for writability if results are being saved in image header
	    if (savebcv) {
		if (imaccess (specpath, READ_WRITE) == NO) {
		call eprintf ("BCVCORR: cannot write to %s; not saving results\n")
		    call pargstr (specpath)
		    IM_UPDATE(specim) = NO
		    }
		else {

#		Write midtime of observation in UT
		    utm = dindef
		    call imgdpar (specim,"UTMID   ", utm)
		    if (utm == dindef) {
			call sprintf (uts, 32, "%h")
			    call pargd (ut)
			call imastr (specim, "UTMID", uts)
			if (debug) {
			    call eprintf ("BCVCORR: writing UTMID = %.1h to %s\n")
				call pargd (uts)
				call pargstr (specpath)
			    }
			}
		    else {
			if (debug) {
			    call eprintf ("BCVCORR: writing UTMID = %.1h to %s\n")
				call pargd (utm)
				call pargstr (specpath)
			    }
			}

#		Write midtime of observation as geocentric Julian Date
		    if (imaccf (specim, "GJDN") == NO) {
			call imaddd (specim,"GJDN",gjd)
			if (debug) {
			    call eprintf ("BCVCORR: writing GJDN = %.5f to %s\n")
				call pargd (gjd)
				call pargstr (specpath)
			    }
			}

#		Write midtime of observation as Heliocentric Julian Date
		    if (imaccf (specim, "HJDN") == NO) {
			call imaddd (specim,"HJDN",hjd)
			if (debug) {
			    call eprintf ("BCVCORR: writing HJDN = %.5f to %s\n")
				call pargd (hjd)
				call pargstr (specpath)
			    }
			}

#		Write barycentric velocity correction
		    if (mspec == 0)
			call imaddr (specim,"BCV     ",rbcv)
		    else {
			call sprintf (keyword,SZ_HKWORD,"APVEL%d")
			    call pargi (mspec)
			call sprintf (str,SZ_LINE," 0.0 0.0 _ %.3f")
			    call pargr (rbcv)
			call imastr (specim,keyword,str)
			}
		    if (debug) {
			call eprintf ("BCVCORR: writing BCV = %.4f to %s\n")
			    call pargd (hjd)
			    call pargstr (specpath)
			}
		    IM_UPDATE(specim) = YES
		    }
		}
	    else if (specim != ERR && specim != NULL)
		IM_UPDATE(specim) = NO
	    }

# Close current image and move on to next image
	if (specim != ERR && specim != NULL)
	    call imunmap (specim)
	if (lspec > 0)
	    go to newspec_

#  Close spectrum list
endvc_	if (lspec > 0)
	    call imtclose (speclist)

end

double procedure juldat (mm,dd,yyyy)

int	mm
int	dd
int	yyyy

double	julday
int	igreg, jy, jm, ja

begin

        igreg = 15 + 31 * (10 + 12 * 1582)

	if (yyyy < 0)
	    return 0.d0

	if (mm > 2) {
	    jy = yyyy
	    jm = mm + 1
	    }
	else {
	    jy = yyyy - 1
	    jm = mm + 13
	    }

	julday = int (365.25*jy) + int (30.6001*jm) + dd + 1720995

	if (dd+31*(mm+12*yyyy) >= igreg) {
            ja = int (0.01 * jy)
            julday = julday + 2 - ja + int (0.25 * ja)
            }
	return julday
end


procedure jd2ut (jd, yyyy, mm, dd, ut)

double	jd	#Julian Date
int	yyyy	# Year	(returned)
int	mm	# Month	(returned)
int	dd	# Day	(returned)
double	ut	# UT	(returned)

double	tsec	# Seconds since 1/1/1950 0:00
int	hr	# Hours
int	mn	# Minutes
double	sec	# Seconds

double	t, days
int	ihms,nc,nc4,nly,ny,m,im
define	setmonth_ 10

begin

# Convert Julian Date to seconds since 1950 Jan 1 0:00
	tsec = (jd - 2433282.5d0) * 8.64d4

#  Calculate time of day (hours, minutes, seconds, .1 msec)
	t = dint ((tsec + 61530883200.d0) * 1.d4 + .5d0)
	hr = idint (dmod (t/36.d6,24.d0))
	mn = idint (dmod (t/60.d4,60.d0))
	if (tsec >= 0) {
	    ihms = idint (dmod (tsec+1.d-6,1.d0) * 1.d4)
	    sec = dmod (tsec+1.d-6,60.d0)
	    }
	else {
	    ihms = idint (dmod (tsec-1.d-6,1.d0) * 1.d4)
	    sec = dmod (tsec-1.d-6,60.d0)
	    }
	ut = double (hr) + (double (mn) / 60.d0) + (sec / 3.6d3)

#  Calculate number of days since 0 hr 0/0/0000
	days = dint ((t / 864.0d6) + 1.0d-6)

#  Number of leap centuries (400 years)
	nc4 = idint ((days / 146097.0d0) + 1.0d-5)

#  Number of centuries since last /400
	days = days - (146097.d0 * double (nc4))
	nc = idint ((days / 36524.d0) + 1.0d-4)
	If (nc > 3) nc = 3

#  Number of leap years since last century
	days = days - (36524.d0 * nc)
	nly = idint ((days / 1461.d0) + 1.0d-10)

#  Number of years since last leap year
	days = days - (1461.d0 * double (nly))
	ny = idint ((days / 365.d0) + 1.d-8)
	If (ny > 3) ny = 3

#  Calculate day of month and month
	days = days - (365.d0 * double (ny))
	dd = idint (days + 1.0d-8) + 1
	mm = 0
	do m = 1, 12 {
	    im = mod ((m + ((m - 1) / 5)), 2)
	    if (dd - 1 < im + 30)
		go to setmonth_
	    dd = dd - im - 30
	    }
setmonth_
	mm  =  mod (m+1, 12) + 1

# Year
	yyyy = nc4*400 + nc*100 + nly*4 + ny + m/11

	return
end

# Jul 15 1995	New task
# Oct 18 1995	If observatory position not indef, use parameters
# Oct 18 1995	Always print RA, Dec, Equinox if verbose
# Dec  4 1995	Initialize altitude to 0 if reading from keyword

# Jan 12 1996	Allow position and time entry from parameter file, not image
# May 28 1996	If SPECNUM > 0, write to APVELnnn instead of BCV keyword

# Jan  8 1997	Compute heliocentric julian date
# Jan 17 1997	Compute BCV at geocentric Julian Date
# Jan 30 1997	Minimize repition in output
# Oct 23 1997	Fix egregious error which reversed BCV and HCV in call to BCV
# Dec 17 1997   Change order of date from dd-mm-yyyy to yyyy-mmm-dd
# Dec 17 1997   Print names of Julian Date keywords, if written to header
# Dec 17 1997	Always write keywords expected by other RVSAO programs

# Feb 11 1998	Separate date and time existence checks
# Feb 11 1998	Use end or start of observation if center not computable
# Feb 18 1998	Print values when verifying writing to header
# Mar  4 1998	Compute UT if only Julian date is given

# Mar 29 1999	Read IRAF observatory database if name keyword but not lat,long
# Dec  8 1999	Fix bug which failed to use image observatory name