include	"vpthit.h"
include	<imhdr.h>
define	SZ_EXTN		3
#define	SZ_PKTTM	24

#  vpt_do -- Remove particle events for high speed photometer files
#
#  Description:
#  ------------
#  
#  Input CL parameters:
#  -----------------
#
#  Date		Author			Description
#  ----		------			-----------
#  29-Oct-1990  J.-C. Hsu		rewrite in SPP
#------------------------------------------------------------------------------

procedure vpt_do (fin, outflag, fout, outdir, scheme, 
			sigmas, dthresh, athresh, lowflag, hitval, 
			root, mean, err, remove, samptm, type, idroot, epoch)

						## input
pointer	fin, fout	 			# file template pointer
int	outflag
char	outdir[SZ_PATHNAME]
char	scheme[SZ_SCHEME]
real	dthresh, athresh, sigmas, hitval
bool	lowflag
						## output
char	root[SZ_ROOT]
real	mean, err
int	remove
real	samptm
int	type
int	idroot
double	epoch
						## local
int	mask
char	ifile[SZ_FNAME], ofile[SZ_FNAME], imask[SZ_FNAME]
pointer	ipin, ipout, ipinmask
pointer	sp, ip
char	typestr[SZ_BUNIT]
pointer	arrdata, arrmask
char	cluster[SZ_FNAME]
char	inroot[SZ_FNAME], outroot[SZ_FNAME]
char	extn[SZ_EXTN]
int	i, j, k
char	dirname1[SZ_PATHNAME]
#char	pkttm[SZ_PKTTM]
int	npix
int	nchar
int	root_len

pointer	immap() 
pointer	imgl1r(), imgs2r(), impl1r(), imps2r()
real	imgetr()
double	imgetd()
int	imtgetim(), fnldir()
int	strsearch()
int	imaccess()
int	ctoi()
#==============================================================================
begin
	    # read the next file name in the template list
	    nchar = imtgetim (fin, ifile, SZ_FNAME)

	    call printf("Processing file %s\n")
	        call pargstr(ifile)

	    if (outflag == NULL)
		ofile[1] = EOS
	    else if (outflag == DIR) {

		# if the output is a directory name, use the root name of the
		# input file and attach the extension of .n*
		# (adapted from t_imcopy.x in images$)
	        call imgcluster (ifile, cluster, SZ_FNAME)
	    	call iki_parse (cluster, inroot, extn)
		root_len = fnldir (inroot, dirname1, SZ_PATHNAME)
		call strcpy (inroot[root_len + 1], dirname1, SZ_PATHNAME)
		call strcpy (outdir, outroot, SZ_PATHNAME)
		call strcat (dirname1, outroot, SZ_PATHNAME)
	    	extn[1] = 'n'
	    	call iki_mkfname (outroot, extn, ofile, SZ_FNAME)

	    } else if (outflag == TEMPL) 
	   	nchar = imtgetim (fout, ofile, SZ_FNAME)

	    # assume the mask file has the same as the input file except the
	    # extension is .q* instead of .d*
	    call imgcluster (ifile, cluster, SZ_FNAME)
	    call iki_parse (cluster, inroot, extn)
	    k = 2
	    j = ctoi (extn, k, idroot)
	    if (idroot < 0 || idroot > 3) {
		call eprintf ("illegal input file extension %s\n")
		    call pargstr (extn)
	    }

	    extn[1] = 'q'
	    call iki_mkfname (inroot, extn, imask, SZ_FNAME)

	    mask = imaccess(imask, READ_ONLY)

	    # open input file and find out how many points are there
	    ipin = immap (ifile, READ_ONLY, 0)
	    npix = 1
	    do j = 1, IM_NDIM(ipin)
	        npix = npix * IM_LEN(ipin, j)

	    # allocate memory for output mask
	    call smark (sp)
	    call salloc (arrmask, npix, TY_REAL)

	    # open input mask, if there is any, and check its size
	    if (mask == YES) {
	        ipinmask = immap (imask, READ_ONLY, 0)
		if (IM_NDIM(ipinmask) != IM_NDIM(ipin)) {
		    call eprintf (
		       	"input mask and file no. %d differ in dimension\n")
			call pargi (i)
		    call error (1, "")
		}
		do j = 1, IM_NDIM(ipinmask) {
		    if (IM_LEN(ipin, j) != IM_LEN(ipinmask, j)) {
		    	call eprintf ( 
				"input mask and file no. %d differ in size\n")
			    call pargi (i)
			call error (1, "")
		    }
		}
	    } else 
		call amovkr (real(OKVAL), Memr[arrmask], npix)

	    # read data from input file
	    if (IM_NDIM(ipin) == 1) {
	        arrdata = imgl1r (ipin)
		if (mask == YES)
	            arrmask = imgl1r (ipinmask)
	    } else if (IM_NDIM(ipin) == 2) {
	        arrdata = imgs2r (ipin, 1, IM_LEN(ipin, 1), 1, IM_LEN(ipin, 2))
		if (mask == YES)
	            arrmask = imgs2r (ipinmask, 1, IM_LEN(ipin, 1), 1, 
					IM_LEN(ipin, 2))
	    } else 
		call error (1, "incorrect dimension (>2)")

	    # read keywords from input file header
	    samptm = imgetr (ipin, "SAMPTIME")

	    #call imgstr (ipin, "FPKTTIME", pkttm, SZ_PKTTM)
	    #call pkttime (pkttm, epoch)
	    # ICD-19 change packet time to MJD
	    epoch = imgetd (ipin, "FPKTTIME")

	    call imgstr (ipin, "ROOTNAME", root, SZ_ROOT)
	    call imgstr (ipin, "DATA_TYP", typestr, SZ_BUNIT)
	    if (strsearch (typestr, "ANALOG") > 0)
		type = ANALOG
	    else if (strsearch (typestr, "DIGITAL") > 0)
		type = DIGITAL
	    else
		call error (1, "unknown data type")

	    # remove particle events
	    call vpt_remove (Memr[arrdata], Memr[arrmask], scheme, npix, 
				type, samptm, sigmas, dthresh, athresh, lowflag,
				hitval, mean, err, remove)

	    # write to output mask file
	    if (outflag != NULL) {
		ipout = immap (ofile, NEW_COPY, ipin)
	        if (IM_NDIM(ipin) == 1)
		    ip = impl1r (ipout)
	        else if (IM_NDIM(ipin) == 2)
	            ip = imps2r (ipout, 1, IM_LEN(ipin, 1), 1, IM_LEN(ipin, 2))
		call amovr (Memr[arrmask], Memr[ip], npix)
		call imunmap (ipout)
	    }

	    # close input file(s)
	    call imunmap (ipin)
	    if (mask == YES)
		call imunmap (ipinmask)

	    call sfree (sp)
end