include <imhdr.h>
include <mach.h>
define	SZ_GEISTIME	24
define 	SEC_PER_DAY	86400.
define	SZ_EPSPEC	3
define	MAX_RATIO	1.d10

#  GET_PHASE -- construct the phase array for a time series observation.
#  First of all, we need to read or construct an epoch array for the time 
#  series observation.  There are two options:  if the calling program is 
#  supplying a time file (e.g. from the task delaytime), then the output 
#  epoch array is simply the supplied time file.  If there is no time file, the 
#  header information of exposure start time (EXPSTART) and CD1_1 are used to 
#  construct a linear epoch array. The phase array is simply the time array 
#  divided by the period.  Only the fractional phases remain in the phase array.
#  
#  input header keywords
#  ---------------------
#  "REFMJD"		reference MJD of the time file
#  "EXPSTART"		exposure start time 
#  "CD1_1"		time span between pixels in the input data file
#
#  Date		Author			Description
#  ----		------			-----------
#  25-Jan-1990  J.-C. Hsu		design and code
#  23-Jun-1994  J.-C. Hsu		change from FPKTTIME to EXPSTART (v1.1)
#------------------------------------------------------------------------------

procedure get_phase (timeflag, timefile, ipin, phase, 
			period, epoch0, epoch0_spec, npix)

bool	timeflag		# input: is there an input time file?
char	timefile[SZ_FNAME]      # input: time data file name
pointer	ipin			# input: file pointer of the input science file
double	phase[npix]		# output: phase array
double	period			# input: period of the light curve
double	epoch0			# input: epoch of zero phase
char	epoch0_spec[SZ_EPSPEC]	# input: specification of epoch0
int	npix			# input: number of points in the time series

char	str[SZ_LINE]
int	pid, i
pointer	iptime
pointer ip			# dummy pointer
	
#char	fpkttime[SZ_GEISTIME]
real	cd1_1			# time between pixels, in seconds
double	mjd0			# modified Julian day of the first packet

pointer	immap()
pointer	imgl1d()
real	imgetr()
double	imgetd()
bool	streq()
#==============================================================================
begin

	# case 1: there is an input time file
	if (timeflag) {

	    # open the input time file, read data from the time file 
	    # time file must be 1-D data 	    
	    # size of the time file must be the same as the input file
	    iptime = immap (timefile, READ_ONLY, 0)		
	    if (IM_NDIM(iptime) != 1) 		
		call error (1, "time file is not 1-D")
	    if (IM_LEN(iptime, 1) != npix) 		
		call error (1, "time file has different size than input file")

            # read the header keywords of REFMJD from the input time file
	    mjd0 = imgetd (iptime, "REFMJD")

	    ip = imgl1d (iptime)
	    call amovd (Memd[ip], phase[1], npix)
	    do i = 1, npix
		phase[i] = mjd0 + phase[i] / SEC_PER_DAY

	    # close the time file
	    call imunmap (iptime)

        # case 2: no input time file
        } else {

            # read the header keywords (group parameters) of CD1_1 and EXPSTART
	    # from the input science file
	    # chabge from PKTTIME to EXPSTART, 6/23/94 JC Hsu.
	    #call imgstr (ipin, "FPKTTIME", fpkttime, SZ_GEISTIME)
	    mjd0 = imgetd (ipin, "EXPSTART")
	    cd1_1 = imgetr (ipin, "CD1_1")

	    #call geistime_mjd (fpkttime, mjd0)	    
	    do i = 1, npix
		phase[i] = mjd0 + double(i-1) * double(cd1_1 / SEC_PER_DAY)
	}

	#  calculate the phases
	str[1] = EOS
	call strlwr (epoch0_spec)

	# case 1: the input epoch of zero phase is MJD
	if (streq(epoch0_spec, "mjd")) 
	    ;

	# case 2: the input epoch of zero phase is pixel number (can be 
	# non-integer) of the input file
	else if (streq(epoch0_spec, "pix")) {
	    pid = int(epoch0)
	    if (epoch0 < 1.d0 || epoch0 > double(npix)+EPSILOND)
		call error (1, "epoch of zero phase is out of pixel range")
	    else if (pid == npix)
		epoch0 = phase[npix]

	    # interpolate between the neighboring time values
	    else {
		epoch0 = phase[pid] + (epoch0 - double(pid)) * (phase[pid+1] -
				phase[pid])
	    }
	} else {
	    call strcat ("illegal epoch specification: ", str, SZ_LINE)
	    call strcat (epoch0_spec, str, SZ_LINE)
	    call error (1, str)
	}

	# if the zero phase epoch is too far from the epoch of the first
	# pixel, issue warning messages
	if ((abs(phase[1]-epoch0)/period) > MAX_RATIO)
	    call printf("Warning: epoch0 is too far from the data set epoch\n")

	# calculate phases, all time quantities are in DAYS
	do i = 1, npix {
	    phase[i] = mod((phase[i] - epoch0), period) / period
	    if (phase[i] < 0.d0)
		phase[i] = phase[i] + 1.d0
	    if (phase[i] == 1.d0)
		phase[i] = 0.d0
	}
end