include <imio.h>
include <imhdr.h>
include <mef.h>
include <ctype.h>

define MAX_RANGES  100

         
# FITSCOPY --  Accumulate 2 or more input FITS file into one multiple
# extension file or append several FITS units to another FITS file

procedure t_fxcopy()

char	input_file[SZ_PATHNAME]
char	output_file[SZ_PATHNAME]
char	imtlist[SZ_LINE], lbrk		

char	group_list[SZ_LINE]
int	group_range[2*MAX_RANGES+1]

int	nch, in, out, gn, ninfiles, nitems, list
int	imtopen(), imtgetim(), imtlen()
int	lget_next_number(), ldecode_ranges(), stridx()
     
pointer mefi, mefo, fco_open_output(), mef_open()
bool    clgetb(), new_file, ogrp
bool    verbose
define  err_ 99
errchk  fcopy

begin

	call clgstr ("input", imtlist, SZ_LINE)
	call clgstr ("output", output_file, SZ_LINE)
	new_file = clgetb("new_file")
	verbose = clgetb("verbose")

	mefo =  fco_open_output (output_file, new_file)

	out = MEF_FD(mefo)

	list = imtopen (imtlist)
	ninfiles = imtlen (list)
	
	lbrk = '['
	if (ninfiles == 1)
	   ogrp = (stridx (lbrk, imtlist) > 0)

	# Accumulate 2 or more input FITS files into one mef file or copy
	# one group from one input file.

	if (ninfiles > 1 || ogrp) {
	    while (imtgetim (list, input_file, SZ_PATHNAME) != EOF) {
		mefi = mef_open(input_file, READ_ONLY, 0)

        	# If verbose print the operation.
		if (verbose) {
		    call eprintf ("%s -> %s\n")
		      call pargstr (input_file)
		      call pargstr (output_file)
		}

        	# Copy input file to either new or existent file.
		if (new_file) {
		   if (MEF_ENUMBER(mefi) >= 0)
		       call mef_copy_extn (mefi, mefo, MEF_ENUMBER(mefi))
	 	   else
		       call fcopyo (MEF_FD(mefi), MEF_FD(mefo))
		   
		   MEF_ACMODE(mefo) = APPEND
		   new_file = false
	        } else {
		   if (MEF_ENUMBER(mefi) >= 0)
		       call mef_copy_extn (mefi, mefo, MEF_ENUMBER(mefi))
		   else
		       call mef_app_file (mefi, mefo)
	        }
		call mef_close (mefi)
	    }
            call mef_close (mefo)
	    call imtclose (list)

	} else {
	    # Copy selected extensions from one input file to the output file.
            call clgstr ("groups", group_list, SZ_LINE)

	    # Since ranges handles only positive numbers, see if we need
	    # to look for zero group (PHDU).

            if (ldecode_ranges (group_list,group_range,MAX_RANGES,nitems) ==ERR)
	        call error (0, "Illegal file number list.")

	    list = imtopen (imtlist)

	    nch = imtgetim (list, input_file, SZ_PATHNAME) 

	    mefi = mef_open (input_file, READ_ONLY, 0)

	    in = MEF_FD(mefi)
	    
	    # If no input group list is specified, copy the whole file.
	    if (group_list[1] == EOS) {
	       if (verbose) {
		  call eprintf ("%s -> %s\n")
		    call pargstr (input_file)
		    call pargstr (output_file)
	       }
	       call fcopyo (in, out)
	    } else {
	       gn = -1
               while (lget_next_number (group_range, gn) != EOF) {
	           call mef_copy_extn (mefi, mefo, gn)
	           if (verbose) {
		      call eprintf ("%s[%d] -> %s\n")
		      call pargstr (input_file)
		      call pargi(gn)
		      call pargstr (output_file)
	           }
	       }
	    }	
err_
	    call close (in)
	    call close (out)

	    call imtclose (list)
	}
end


# FCOP_OPEN_OUTPUT -- Open output file and return the mef descriptor.

pointer procedure fco_open_output (output_file, new_file)

char	output_file[ARB]	#I, output filename
bool	new_file			#I, true if file already exists

pointer mef, mef_open()
int	access(), acmode
errchk  mef_open

begin
	if (!new_file) {
	   # See if the file exists else change mode.
	   if (access (output_file, 0, 0) == NO)
	      new_file = true
	}
	
	acmode = APPEND
	if (new_file) {
           call fclobber (output_file)
	   acmode = NEW_FILE
	}

	mef = mef_open (output_file, acmode, 0)

	if (MEF_ENUMBER(mef) != -1)
	   call error(13, "Extension number not allowed in filename")

	return (mef)
end


# NOTE: This local version of the xtools routine call handle starting
# index of zero (0). Taken from dataio/lib and modified. NZ March, 98
#
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<mach.h>
include	<ctype.h>

define	FIRST	1		# Default starting range
define	LAST	MAX_INT		# Default ending range
define	STEP	1		# Default step
define	NULL	-1		# Ranges delimiter

# DECODE_RANGES -- Parse a string containing a list of integer numbers or
# ranges, delimited by either spaces or commas.  Return as output a list
# of ranges defining a list of numbers, and the count of list numbers.
# Range limits must be positive nonnegative integers.  ERR is returned as
# the function value if a conversion error occurs.  The list of ranges is
# delimited by a single NULL.

int procedure ldecode_ranges (range_string, ranges, max_ranges, nvalues)

char	range_string[ARB]	# Range string to be decoded
int	ranges[3, max_ranges]	# Range array
int	max_ranges		# Maximum number of ranges
int	nvalues			# The number of values in the ranges

int	ip, nrange, first, last, step, ctoi()

begin
	ip = 1
	nvalues = 0

	do nrange = 1, max_ranges - 1 {
	    # Defaults to all positive integers
	    first = FIRST
	    last = LAST
	    step = STEP

	    # Skip delimiters
	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
		ip = ip + 1

	    # Get first limit.
	    # Must be a number, '-', 'x', or EOS.  If not return ERR.
	    if (range_string[ip] == EOS) {			# end of list
		if (nrange == 1) {
		    # Null string defaults
		    ranges[1, 1] = first
		    ranges[2, 1] = last
		    ranges[3, 1] = step
		    ranges[1, 2] = NULL
	    	    nvalues = nvalues + abs (last-first) / step + 1
		    return (OK)
		} else {
		    ranges[1, nrange] = NULL
		    return (OK)
		}
	    } else if (range_string[ip] == '-')
		;
	    else if (range_string[ip] == 'x')
		;
	    else if (IS_DIGIT(range_string[ip])) {		# ,n..
		if (ctoi (range_string, ip, first) == 0)
		    return (ERR)
	    } else
		return (ERR)

	    # Skip delimiters
	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
		ip = ip + 1

	    # Get last limit
	    # Must be '-', or 'x' otherwise last = first.
	    if (range_string[ip] == 'x')
		;
	    else if (range_string[ip] == '-') {
		ip = ip + 1
	        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
		    ip = ip + 1
		if (range_string[ip] == EOS)
		    ;
		else if (IS_DIGIT(range_string[ip])) {
		    if (ctoi (range_string, ip, last) == 0)
		        return (ERR)
		} else if (range_string[ip] == 'x')
		    ;
		else
		    return (ERR)
	    } else
		last = first

	    # Skip delimiters
	    while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
		ip = ip + 1

	    # Get step.
	    # Must be 'x' or assume default step.
	    if (range_string[ip] == 'x') {
		ip = ip + 1
	        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
		    ip = ip + 1
		if (range_string[ip] == EOS)
		    ;
		else if (IS_DIGIT(range_string[ip])) {
		    if (ctoi (range_string, ip, step) == 0)
		        ;
		} else if (range_string[ip] == '-')
		    ;
		else
		    return (ERR)
	    }

	    # Output the range triple.
	    ranges[1, nrange] = first
	    ranges[2, nrange] = last
	    ranges[3, nrange] = step
	    nvalues = nvalues + abs (last-first) / step + 1
	}

	return (ERR)					# ran out of space
end


# GET_NEXT_NUMBER -- Given a list of ranges and the current file number,
# find and return the next file number.  Selection is done in such a way
# that list numbers are always returned in monotonically increasing order,
# regardless of the order in which the ranges are given.  Duplicate entries
# are ignored.  EOF is returned at the end of the list.

int procedure lget_next_number (ranges, number)

int	ranges[ARB]		# Range array
int	number			# Both input and output parameter

int	ip, first, last, step, next_number, remainder

begin
	# If number+1 is anywhere in the list, that is the next number,
	# otherwise the next number is the smallest number in the list which
	# is greater than number+1.

	number = number + 1
	next_number = MAX_INT

	for (ip=1;  ranges[ip] != NULL;  ip=ip+3) {
	    first = min (ranges[ip], ranges[ip+1])
	    last = max (ranges[ip], ranges[ip+1])
	    step = ranges[ip+2]
	    if (number >= first && number <= last) {
		remainder = mod (number - first, step)
		if (remainder == 0)
		    return (number)
		if (number - remainder + step <= last)
		    next_number = number - remainder + step
	    } else if (first > number)
		next_number = min (next_number, first)
	}

	if (next_number == MAX_INT)
	    return (EOF)
	else {
	    number = next_number
	    return (number)
	}
end


# GET_PREVIOUS_NUMBER -- Given a list of ranges and the current file number,
# find and return the previous file number.  Selection is done in such a way
# that list numbers are always returned in monotonically decreasing order,
# regardless of the order in which the ranges are given.  Duplicate entries
# are ignored.  EOF is returned at the end of the list.

int procedure lget_previous_number (ranges, number)

int	ranges[ARB]		# Range array
int	number			# Both input and output parameter

int	ip, first, last, step, next_number, remainder

begin
	# If number-1 is anywhere in the list, that is the previous number,
	# otherwise the previous number is the largest number in the list which
	# is less than number-1.

	number = number - 1
	next_number = 0

	for (ip=1;  ranges[ip] != NULL;  ip=ip+3) {
	    first = min (ranges[ip], ranges[ip+1])
	    last = max (ranges[ip], ranges[ip+1])
	    step = ranges[ip+2]
	    if (number >= first && number <= last) {
		remainder = mod (number - first, step)
		if (remainder == 0)
		    return (number)
		if (number - remainder >= first)
		    next_number = number - remainder
	    } else if (last < number) {
		remainder = mod (last - first, step)
		if (remainder == 0)
		    next_number = max (next_number, last)
		else if (last - remainder >= first)
		    next_number = max (next_number, last - remainder)
	    }
	}

	if (next_number == 0)
	    return (EOF)
	else {
	    number = next_number
	    return (number)
	}
end


# IS_IN_RANGE -- Test number to see if it is in range.

bool procedure lis_in_range (ranges, number)

int	ranges[ARB]		# Range array
int	number			# Number to be tested against ranges

int	ip, first, last, step

begin
	for (ip=1;  ranges[ip] != NULL;  ip=ip+3) {
	    first = min (ranges[ip], ranges[ip+1])
	    last = max (ranges[ip], ranges[ip+1])
	    step = ranges[ip+2]
	    if (number >= first && number <= last)
		if (mod (number - first, step) == 0)
		    return (true)
	}

	return (false)
end