include <imhdr.h>
include	<time.h>
define	SZ_SPEC		3

#  ORDERPH_NTO1 -- Order phase from N input files to one output file.
#
#  Description:
#  ------------
#
#  Date		Author			Description
#  ----		------			-----------
#  25-Jan-1990  J.-C. Hsu		rewrite in SPP
#------------------------------------------------------------------------------

procedure orderph_nto1 (fin, fout, fphase, ftime, finmask, foutmask,
			nfin, period, epoch0, epoch0_spec, maskflag, timeflag)

pointer	fin, fout, fphase, ftime, finmask, foutmask	
				# input: file template pointers
int	nfin 			# input: number of files in the input template
double	period			# input: period of the light curve
double	epoch0			# input: epoch of the zero phase
char	epoch0_spec[SZ_SPEC]	# how is zero phase epoch specified?
bool	maskflag		# input: is there input mask file(s)?
bool	timeflag		# input: is there input time file(s)?

pointer arrdata, arrphase, arrmask
				# addresses of data
pointer	ptdata, ptphase, ptmask
  
pointer	ipin, ipout, ipinmask, ipoutmask, ipphase 
char	ifile[SZ_FNAME], ofile[SZ_FNAME], imask[SZ_FNAME], omask[SZ_FNAME],
	phfile[SZ_FNAME], timefile[SZ_FNAME]
int	npix			# number of points in the time series
int	npts
int	nchar, i
int	hsize
pointer sp, ip			# dummy pointer
pointer	text
	
pointer	immap()
pointer	imgl1r(), impl1r(), impl1d() 
int	imtgetim()
long	clktime()
#==============================================================================
begin
        call smark (sp)

	# find out how long the history log may be
	hsize = SZ_TIME + (2 * SZ_FNAME + 30) * nfin + 50
	call salloc (text, hsize, TY_CHAR)

	# history message
	call cnvtime (clktime(0), Memc[text], hsize)
	call strcat (" created by the task ORDERPHASE, from: ", 
			Memc[text], hsize)

        nchar = imtgetim (fout, ofile, SZ_FNAME) 
        nchar = imtgetim (fphase, phfile, SZ_FNAME) 
	if (maskflag) 
            nchar = imtgetim (foutmask, omask, SZ_FNAME) 

	npts = 0
	# find out how many points are in all the files and allocate memory
	do i = 1, nfin {	
	    nchar = imtgetim (fin, ifile, SZ_FNAME) 
	    ipin = immap (ifile, READ_ONLY, 0)		
	    npts = IM_LEN(ipin, 1) + npts
	    call imunmap (ipin)
	}
	call salloc (arrphase, npts, TY_DOUBLE)
	call salloc (arrdata, npts, TY_REAL)
	ptdata = arrdata
	ptphase = arrphase
	if (maskflag) { 
	    call salloc (arrmask, npts, TY_REAL)
	    ptmask = arrmask
 	}

	# rewind the template pointer and loop all input files
	call imtrew (fin)
	do i = 1, nfin {
	    # read the next file name in the template list
	    nchar = imtgetim (fin, ifile, SZ_FNAME) 
	    if (timeflag)
		nchar = imtgetim (ftime, timefile, SZ_FNAME)

	    # open input files
	    # input file must be 1-D data and real data type
	    ipin = immap (ifile, READ_ONLY, 0)		
	    if (IM_NDIM(ipin) != 1 || IM_PIXTYPE(ipin) != TY_REAL) 	
		call error (1, "input is not 1-D data or is not real data type")

	    # open output files
	    if (i == 1) {
	    	ipout = immap (ofile, NEW_COPY, ipin)	
	    	ipphase = immap (phfile, NEW_COPY, ipin)	
		IM_PIXTYPE(ipphase) = TY_DOUBLE
		IM_LEN(ipphase, 1) = npts
		IM_LEN(ipout, 1) = npts

		# put a dummy element to the output files so we don't have 
		# to keep the input file open till the final write
		Memd[arrphase] = 0.
		ip = impl1d (ipphase)
		call amovd (Memd[arrphase], Memd[ip], 1)
		Memr[arrdata] = 0.
		ip = impl1r (ipout)
		call amovr (Memr[arrdata], Memr[ip], 1)
	    }

	    # read data from the input file 
	    npix = IM_LEN(ipin, 1)
	    ip = imgl1r (ipin)
	    call amovr (Memr[ip], Memr[ptdata], npix)

	    # if there is input mask, read its data
	    if (maskflag) {
	        nchar = imtgetim (finmask, imask, SZ_FNAME) 
	    	ipinmask = immap (imask, READ_ONLY, 0)
		if (IM_LEN(ipinmask, 1) != npix) 
		    call error (1, 
			  "mismatch of file size between input and input mask")
		# open output mask file
		if (i == 1) {
	    	    ipoutmask = immap (omask, NEW_COPY, ipin)
		    IM_LEN(ipoutmask, 1) = npts
		    # put a dummy element to the output files so we don't 
		    # have to keep the input file open till the final write
		    Memr[arrmask] = 0.
		    ip = impl1r (ipoutmask)
		    call amovr (Memr[arrmask], Memr[ip], 1)
		}
	    	ip = imgl1r (ipinmask)
	    	call amovr (Memr[ip], Memr[ptmask], npix)
		ptmask = ptmask + npix
		call imunmap (ipinmask)
	    }
	    
	    # get the phase array 
	    # if the epoch0 specification is pix, it is only for the first file
	    if (i > 1)
		call strcpy ("mjd", epoch0_spec, SZ_SPEC)
	    call get_phase (timeflag, timefile, ipin, Memd[ptphase], 
				period, epoch0, epoch0_spec, npix)

	    ptdata = ptdata + npix
	    ptphase = ptphase + npix
	    call imunmap (ipin)
	    call strcat (" ", Memc[text], hsize)
	    call strcat (ifile, Memc[text], hsize)
	    if (timeflag) {
		call strcat (" and its time file ", Memc[text], hsize)
		call strcat (timefile, Memc[text], hsize)
		call strcat (";", Memc[text], hsize)
	    }
	}

	# Sort phases, data, and mask values
	call sort_phase (Memd[arrphase], Memr[arrdata], Memr[arrmask], 
					maskflag, npts)

	# put the result to the output files 
	ip = impl1r (ipout)
	call amovr (Memr[arrdata], Memr[ip], npts)
	ip = impl1d (ipphase)
	call amovd (Memd[arrphase], Memd[ip], npts)
	if (maskflag) {
	    ip = impl1r (ipoutmask)
	    call amovr (Memr[arrmask], Memr[ip], npts)

	    call imunmap (ipoutmask)
	}

	# update header keyword(s) and history
	call impstr (ipphase, "BUNIT", "PHASE", SZ_BUNIT)
	call imputh (ipphase, "HISTORY", Memc[text])
	call imputh (ipout, "HISTORY", Memc[text])

	# close output files
	call imunmap (ipout)
	call imunmap (ipphase)

	# print out message of which files been created
	call printf ("orderphase: output file %s is created\n")
	    call pargstr (ofile)
	call printf ("orderphase: phase file %s is created\n")
	    call pargstr (phfile)
	if (maskflag) {
	    call printf ("orderphase: output mask %s is created\n")
	    	call pargstr (omask)
	}

	call sfree (sp)
end