include	<mach.h>
include	<math.h>
include	<tbset.h>
include	<imhdr.h>

define	NUM_OF_SETS	Memi[$1]
define	IND_OF_SET	Memi[$1+2*$2-1]
define	NUM_IN_SET	Memi[$1+2*$2]

# T_SLITANGLEX -- Calculate angle of spectrograph lines with axis 1 direction
# The program calculates the angle the lines from the long slit spectrograph 
# make with the pixel direction. The image is assumed to have been corrected 
# for PDA distortion. The data consists of a spectrograph image together with 
# a table in which a one or more lines have been delineated by marking points 
# above and below the spectral line at two or more positions across the image.
#
# D. Giaretta	Dec 1987	Original SPP version
# R.N.Hook      Jul 1988        Bug-fix for two points per line case
# Phil Hodge	8-Aug-1993	Swap arguments in call to tbcinf.

procedure t_slitanglex()

char	image[SZ_FNAME]			# input spectrog. image
char	table[SZ_FNAME]			# table name
char	xcol[SZ_COLNAME]		# x col name
char	ycol[SZ_COLNAME]		# y col name
char	setcol[SZ_COLNAME]		# set col name

pointer	im				# image structure
pointer	xpt				# pointer to xpt data
pointer	ypt				# pointer to ypt data
pointer	setpt				# pointer to set data
pointer	isetpt				# pointer groups structure
pointer	xnull, ynull, snull		# null flag pointers
pointer	v1pt, v2pt			# pointers to temporary accumulations
					# of x, y data in one set
pointer	pangle, posx, posy
int	num_used			# number of points used in set
int	j				# Loop index
int	setnum				# Number of data points in current set
int	numrows				# Number fo rows in input file
int	nsum				# 
real	angle				# Final angle
real	sdangle				# SD of angle
double	sum, sum2			# sums for statistics

pointer	immap()
errchk	immap
begin

	# get names of files
	call clgstr( "input", image, SZ_FNAME)
	im = immap( image, READ_ONLY, 0)
	call clgstr( "table", table, SZ_FNAME)
	call clgstr( "xcol"  , xcol  , SZ_COLNAME)
	call clgstr( "ycol"  , ycol  , SZ_COLNAME)
	call clgstr( "setcol", setcol, SZ_COLNAME)

	# read in all the data
	call sl_get_data( table, xcol, ycol, setcol, 
			  numrows, xpt, ypt, setpt, xnull, ynull, snull)

	# allocate space for results
	call salloc( isetpt, 2*numrows+1, TY_INT)
	call salloc( v1pt,  numrows, TY_REAL)
	call salloc( v2pt,  numrows, TY_REAL)
	call salloc( pangle,numrows, TY_REAL)
	call salloc( posx,  numrows, TY_REAL)
	call salloc( posy,  numrows, TY_REAL)

	# find the array of different set numbers
	call sl_set_numbers( Memi[setpt], Memb[snull], numrows, isetpt )

	# loop through sets
	for ( setnum = 1; setnum <= NUM_OF_SETS(isetpt); setnum = setnum + 1) {

	    # read all the data for that set
	    call sl_get_set( Memr[xpt], Memr[ypt], Memi[setpt], 
			Memb[xnull], Memb[ynull], Memb[snull], numrows, 
			setnum, isetpt,
			Memr[v1pt], Memr[v2pt], num_used )

	    call printf( " Data for line %d\n")
	    	call pargi( setnum)
	    call flush (STDOUT)

	    # process data for this set
	    call slit_angcl( im, Memr[v1pt], Memr[v2pt], 
				num_used, 
				Memr[posx], Memr[posy], Memr[pangle+setnum-1])

	}          

	# Calculate the mean angle, and sd - there must be at least one point
	nsum = 0 ; sum  = 0.0 ; sum2 = 0.0
        do j = 0, NUM_OF_SETS(isetpt)-1 {
	    if ( !IS_INDEFR (Memr[pangle+j]) ) {
		nsum = nsum + 1
           	sum  = sum  + Memr[pangle+j]
           	sum2 = sum2 + Memr[pangle+j]**2
	    }
	}

	if ( nsum < 1 )
            call error( 0, "there must be at least one line ")
	else {
            angle = sum/float(nsum)
            if ( nsum > 1) 
        	sdangle = sqrt(sum2 - nsum*angle*angle)/(nsum-1)
            else
           	sdangle = 0.0
        }

	call printf( " Mean angle = %14.7g, SD= %14.7g \n")
	    call pargr( angle) ; call pargr( sdangle)

	# Insert into process parameters
	call clputr( "angle", angle)
	call clputr( "sdangle", sdangle)

	call imunmap(im)

end
# SL_GET_DATA -- read data from table
procedure sl_get_data( table, xcol, ycol, setcol,
			numrows, xpt, ypt, setpt, xnull, ynull, snull)

char	table[ARB]	# i: table name
char	xcol[ARB]	# i: x col name
char	ycol[ARB]	# i: y col name
char	setcol[ARB]	# i: set col name
int	numrows		# o: number of rows
pointer	xpt		# o: x data pointer
pointer	ypt		# o: y data pointer
pointer	setpt		# o: set data pointer
pointer	xnull		# o: x null flag pointer
pointer	ynull		# o: y null flag pointer
pointer	snull		# o: set null flag pointer
#--

pointer	tp, tbtopn()
pointer	xcolpt, ycolpt, setcolpt	# table col pointers
pointer	setchar
char	junkchar[SZ_LINE]
int	junkint
int	datatype, lendata, i
int	tbpsta()
errchk	tbtopn

begin
	tp = tbtopn( table, READ_ONLY, 0)
	numrows	 = tbpsta( tp, TBL_NROWS)

	call salloc( xpt,   numrows, TY_REAL)
	call salloc( xnull, numrows, TY_BOOL)
	call salloc( ypt,   numrows, TY_REAL)
	call salloc( ynull, numrows, TY_BOOL)
	call salloc( setpt, numrows, TY_INT)
	call salloc( snull, numrows, TY_BOOL)

	# read in table data
	call tbcfnd( tp, xcol,   xcolpt, 1)
	call tbcfnd( tp, ycol,   ycolpt, 1)
	call tbcfnd( tp, setcol, setcolpt, 1)

	if (xcolpt == NULL || ycolpt == NULL || setcolpt == NULL)
	    call error (1, "column not found")
	
	# check type of set column
	# Second and fifth arguments interchanged by PEH on 8/13/93.
	call tbcinf( setcolpt, junkint, junkchar, junkchar, junkchar,
			datatype, lendata, junkint)

	# Read data
	call tbcgtr( tp, xcolpt,   Memr[xpt],   Memb[xnull], 1, numrows)
	call tbcgtr( tp, ycolpt,   Memr[ypt],   Memb[ynull], 1, numrows)
	if ( datatype > 0 ) # negative datatype shows character data
	    call tbcgti( tp, setcolpt, Memi[setpt], Memb[snull], 1, numrows)
	else {
	    # if character data we must convert to integer
	    call salloc( setchar, numrows*(lendata+1), TY_CHAR)
	    call tbcgtt( tp, setcolpt, Memc[setchar], Memb[snull], lendata, 
								  1, numrows)
	    do i = 0, numrows-1
		Memi[setpt+i] = Memc[setchar + i*(lendata+1)]
	}
	
	call tbtclo(tp)

end

# SLIT_ANGCL -- Calculate angle of one spectral line
# The program finds the maximum in the spectrum between each pair
# of points, and calculates least squares line through data 
# to define slope of spectrum w.r.t. pixel axis (axis 1).

procedure slit_angcl( im, x, y, npts, posx, posy, pangle)

pointer	im			# i: image structure
real	x[ARB], y[ARB]		# i: x, y data for given set
int	npts			# i: number of points in set
real	posx[ARB], posy[ARB]	# o: positions found, x and y
real	pangle			# o: angle of this spectral line
#--

pointer	bufpt				# image section pointer
double	sum, sumx, sumy, sumxy, sumxx 	# Temp sums
double	sumix, sumiy, sumi		# Temp sums
real	temp
int	ngood				# Number of good values in column
int	ipos				# Number of good positions in 
					#    this line
int	losamp, hisamp, loline, hiline 	# Positions given in table
int	dsamp
int	i, j, k				# Loop index

pointer	imgs2r()
errchk	imgs2r

begin

        ipos = 0

	# Check there are at least 2 points
	if ( npts < 4 ) {
	    pangle = INDEFR
	    call printf( " need at least 2 centroid positions\n")
	    call flush (STDOUT)
	    return
	}

	# Calculate the centroid of the data between the points

	do i = 1, npts, 2 {
	    hisamp = max( x[i], x[i+1] )
	    losamp = min( x[i], x[i+1] )
	    hiline = max( y[i], y[i+1] )
	    loline = min( y[i], y[i+1] )
	    sumix  = 0.0
	    sumiy  = 0.0
	    sumi   = 0.0
	    ngood  = 0

	    call printf ("Image columns :%d to %d ,  lines %d to %d \n")
		call pargi( losamp ); call pargi( hisamp)
		call pargi( loline ); call pargi( hiline)
	    # Check the data lies in image

	    if (losamp >= 1 && losamp <= IM_LEN(im, 1) &&
		hisamp >= 1 && hisamp <= IM_LEN(im, 1) &&
     		loline >= 1 && loline <= IM_LEN(im, 2) &&
		hiline >= 1 && hiline <= IM_LEN(im, 2)    ) {

		dsamp = hisamp - losamp + 1
		bufpt = imgs2r( im, losamp, hisamp, loline, hiline )
		do j = 0, hiline-loline {
		    do k = 0, hisamp-losamp {
			temp = Memr[ bufpt+j*dsamp+k]
		        if (!IS_INDEFR (temp)) {
			    ngood = ngood + 1
			    sumi  = sumi   + temp
			    sumix = sumix + (k+1)*temp
			    sumiy = sumiy + (j+1)*temp
			}
		    }
		}


	        if (ngood == 0 || abs(sumi) < EPSILONR ) {
		    call printf( "no good data points\n")
	        } else {
		    ipos 		= ipos + 1
                    posx[ipos] 	= losamp + sumix/sumi - 1
                    posy[ipos]	= loline + sumiy/sumi - 1
		    call printf (
		       " Centroid of data in column: X = %14.7g, Y = %14.7g \n")
		        call pargr( posx[ipos] ); call pargr( posy[ipos] )
		}
	    } else {

		call printf( "   lies outside image area\n")
	    }     # Are we asking for data inside image?

	}

	# Calculate the mean angle for this spectral line

        if ( ipos < 2 ) {
	    pangle = INDEFR

	} else if ( ipos > 2 ) {
	    sum   = 0.0 ; sumx  = 0.0 ; sumy  = 0.0
	    sumxy = 0.0 ; sumxx = 0.0
	    do i = 1, ipos {
		sum   = sum + 1
		sumx  = sumx + posx[i]
		sumy  = sumy + posy[i]
		sumxy = sumxy + posx[i]*posy[i]
		sumxx = sumxx + posx[i]*posx[i]
	    }

	    pangle = atan2( sumxy*sum - sumx*sumy, sumxx*sum - sumx*sumx)

	} else if ( ipos == 2) {
	
	    # Ensure we are measuring wrt increasing x value
            # Modified to use del Y / del X rather than the other way up
            # Richard Hook, Jul 88
	    if (posx[1] > posx[2] ) 
		pangle = atan2( posy[1] - posy[2], posx[1] - posx[2])
	    else
		pangle = atan2( posy[2] - posy[1], posx[2] - posx[1])
	}

	# convert to degrees
	if ( ipos >= 2)
	    pangle = RADTODEG(pangle)

	call printf( " angle found: %g\n" )
	    call pargr( pangle)

	call flush (STDOUT)

end

# SL_SET_NUMBERS -- find the array of different set numbers

procedure sl_set_numbers( set, setnull, numrows, isetpt )

int	set[ARB]		# i: set index values
bool	setnull[ARB]		# i: are set indices indef.
int	numrows			# i: number of rows
pointer	isetpt			# i: set structure pointers
#--

pointer	tset
int	i, used, numset, setnum, test

begin
	call salloc( tset, numrows, TY_INT)
	used = 0
	for (i=1; i<=numrows; i=i+1) 
	    if ( !setnull[i] ) {
		Memi[tset+used] = set[i]
		used = used + 1
	    }

	call asrti( memi[tset], Memi[tset], used)

	if ( used < 1 )
	    call error( 0, " there must be some non-null data in table")

	setnum	= 1
	numset 	= 0
	test    = Memi[tset]
	for ( i=0; i<=used-1 ; i=i+1) {
	    if ( Memi[tset+i] != test ) {
		IND_OF_SET( isetpt, setnum )  = test
		NUM_IN_SET( isetpt, setnum )  = numset
		numset = 1
		setnum = setnum + 1 
		test   = Memi[tset+i]
	    } else {
		numset = numset + 1
	    }
	}
	IND_OF_SET( isetpt, setnum) = test
	NUM_IN_SET( isetpt, setnum) = numset
	NUM_OF_SETS( isetpt)        = setnum

end

# SL_GET_SET -- read all non-null data fof a given set		

procedure sl_get_set( x, y, set, xn, yn, sn, numrows, setnum, isetpt, 
			v1, v2, num_used )

real	x[ARB]		# i: x data
real	y[ARB]		# i: y data
int	set[ARB]	# i: set data
bool	xn[ARB], yn[ARB], sn[ARB] 	# i: nullflags
int	numrows		# i: number of rows
int	setnum		# i: set number we want
pointer	isetpt		# i: set structure
real	v1[ARB]		# o: x values for this set
real	v2[ARB]		# o: y values for this set
int	num_used	# o: number of points used
#--

int	i, used

begin

	used = 0
	for( i=1; i<=numrows && used < NUM_IN_SET(isetpt, setnum); i=i+1) {
	    if ( !sn[i] && !xn[i] && !yn[i]  && 
		  set[i] == IND_OF_SET(isetpt, setnum) ) {
		used 	 = used + 1
		v1[used] = x[i]
		v2[used] = y[i]
	    }
	}

	num_used = used
end