c ppsfilt.for
c*******************************************************************************
c$ Component_name:
c	PPSFILT
c$ Abstract:
c	This program generates an ASCII-format file containing a filtered,
c	calibrated version of a PPS data file, converted to a radius scale.
c$ Keywords:
c	PROFILE, OCCULTATIONS, PPS
c	FORTRAN, PUBLIC, MAIN_PROGRAM
c$ Declarations:
c	none
c$ Inputs:
c	(1) Name of the PDS label for the PPS data file.
c	(2) Name of the PDS label for the PPS geometry file.
c	(3) Name of the PDS label for the PPS calibration file.
c	(4) Radius curve segment index (if the radius curve is non-monotonic).
c	    Segments are numbered such that segment #1 falls between the first
c	    pair of extrema in the radius function.
c	(5) Lower and upper radius limits (km).
c	(6) Data resampling factor.  For example, a value of 10 will produce a
c	    file with 10 times fewer samples.
c	(7) Point spread function type.  The point spread function (PSF) of the
c	    filtered data can be boxcar, triangle or sinc.  Using a boxcar PSF,
c	    each new sample is a uniform average of the data points centered on
c	    the sample's midpoint and extending halfway to the nearest
c	    neighbors.  Using a triangle PSF, each sample is weighted most
c	    heavily toward the center of the sample and extends (with linearly
c	    decreasing weight) to each of the nearest neighbors.  Using a sinc
c	    PSF, each new sample is best construed as a local sample of a
c	    continuous function that has been low-pass filtered, so that it
c	    contains now power in spatial frequencies higher than the Nyquist
c	    frequency.  Note that the program runs fastest for a boxcar PSF and
c	    perhaps 10 times slower for a sinc PSF.
c	(8) Output file name.
c$ Outputs:
c	A smooth ring profile.  The output file contains six columns:
c	(1) Radius value (km).
c	(2) Mean photon counts.
c	(3) Standard deviation in photon counts.
c	(4) Approximate tau.
c	(5) Lower limit of the 68% (+/-1 sigma) confidence interval on tau.
c	(6) Upper limit of the 68% (+/-1 sigma) confidence interval on tau.
c$ Returns:
c	none
c$ Detailed_description:
c	This program generates an ASCII-format file containing a filtered,
c	calibrated version of a PPS data file, converted to a radius scale.
c	The resultant file can be plotted or used for further analysis.  This
c	program also serves as an example of how to use the Profile Toolkit.
c
c	This program contains a few lines that are non-standard ANSI
c	FORTRAN-77, specifically "IMPLICIT NONE" and the "CARRIAGECONTROL"
c	option on OPEN.  If necessary, these lines can be commented out.
c	However, this appears to be unnecessary on most compilers.
c$ External_references:
c	Profile library, Object Access Library
c$ Examples:
c	Here is a sample run of the program.  It generates output file
c	FILTOUT.TAB, which is included with the toolkit.
c
c	Data file label: PPSDATA.LBL
c	Geometry file label: PPSGEOM.LBL
c	Calibration file label: PPSCALIB.LBL
c	 0: r =  63585.550, dr =    0.109
c	 1: r = 142592.353, dr =    0.087
c	Radius min, max = 133000,134000
c	Resampling factor = 20
c	Filter type (boxcar,triangle,sinc): (B/T/S) = B
c	Output filename: FILTOUT.TAB
c
c	The first five lines of the output file are as follows:
c	133000.00, 23.100, 1.075, 0.4165, 0.3860, 0.4491
c	133001.89, 24.450, 1.106, 0.3785, 0.3494, 0.4094
c	133003.77, 20.450, 1.011, 0.5012, 0.4671, 0.5379
c	133005.66, 20.550, 1.014, 0.4977, 0.4637, 0.5342
c	133007.55, 23.100, 1.075, 0.4164, 0.3859, 0.4490
c$ Error_handling:
c	The program prints an error message and aborts if anything goes wrong.
c$ Limitations:
c	none
c$ Author_and_institution:
c	Mark R. Showalter
c	PDS Rings Node, NASA/Ames Research Center
c$ Version_and_date:
c	1.0: May 1998
c	1.1: December 1999
c	1.2: September 2002
c	1.3: September 2003
c$ Change_history:
c	1.1: Revised to include commas between columns of output file.
c	     Revised to write unweighted samples correctly in output file.
c	1.2: Corrected a bug in the definition of the width of a sinc psf.
c	1.3: Revised so the INCIDENCE_ANGLE comes from the geometry file, not
c	     the calibration file. Updated comments to reflect new directory
c	     structure.
c*******************************************************************************

	program PPSFILT
	implicit none
	include 'fprofile.inc'

	real*8		x, x1, x2, dx, xsize, y, r, r1, r2, drdx,
     &			xr1, xr2, kr1, kr2, center, count, covar, sigma,
     &			b, s, mu, tau, tau_lo, tau_hi
	integer*4	ncolumn, nsegment, nsegments, type, expand,
     &			f1, f2, fsize, k, k1, k2, dk, dksign, flag
	character	filename*80, char*1
	integer*4	dlabel, label, column, series, curve, curve2,
     &			star, back, newpsf, oldpsf, filtered
	real*8		TauFunc

c Constants
	integer		OUTUNIT/1/
	real*8		PI
	real*8		MAXTAU/99.d0/
	real*8		MINTAU/-9.d0/
	real*8		MAXSIGMA/99.d0/

c Filtering parameters
	integer		NFILTER
	parameter	(NFILTER = 500)

	real*8		filter(NFILTER)
	real*8		SINC_HALFWIDTH/5./

c***********************************************************************
c Prepare to set up data series
c***********************************************************************

c Prompt for label filename
	write(*,10) 'Data file label: '
	read(*,11,end=9999) filename
10	format(1x,a,$)
11	format(a)

c Open file
	dlabel = FPro_OpenLabel(filename)

c Create column object (which we'll never use because we're unsure of
c the proper buffer size)
c ntable=1, ncolumn=1, nitem=0 (all), row1=1, row2=0 (all), k1=0,
c buffersize=1, usedouble=FALSE
	series = FPro_ColumnSeries(dlabel, 1, 1, 0, 1, 0, 0, 1, .FALSE.)

c Get interval and prepare to track domain overlap
	xsize = FPro_SeriesSampling(series, x1, x2, dx)
	xsize = FPro_ObjectOverlap(series, series, x1, x2)

c***********************************************************************
c Set up radius function
c***********************************************************************

c Prompt for label filename
	write(*,10) 'Geometry file label: '
	read(*,11,end=9999) filename

c Open file
	label = FPro_OpenLabel(filename)

c Find column
	ncolumn = FPro_LabelFind(label, 1, 'RING_INTERCEPT_RADIUS')

c Create column object
c ntable=1, ncolumn, nitem=1, row1=1, row2=0 (all), k1=0,
c buffersize=0 (all), usedouble=TRUE
	column = FPro_ColumnSeries(label, 1, ncolumn, 1, 1, 0, 0, 0,
     &			.TRUE.)

c Fit cubic splines
	curve = FPro_SplineCurve(column)

c Get mu factor
	PI = 4.d0 * atan(1.d0)
	mu = cos(PI/180.d0 *
     &		FPro_LabelFloat(label, 0, 0, 'INCIDENCE_ANGLE',
     &		0.d0, .TRUE.))

C Free unneeded objects
	call FPro_FreeObject(label)
	call FPro_FreeObject(column)

c Update domain
	xsize = FPro_ObjectOverlap(curve, series, x1, x2)

c***********************************************************************
c Set up stellar and background count functions
c***********************************************************************

c Prompt for label filename
	write(*,10) 'Calibration file label: '
	read(*,11,end=9999) filename

c Open file
	label = FPro_OpenLabel(filename)

c Set up linear splines for stellar signal
	ncolumn = FPro_LabelFind(label, 1, 'STELLAR_SIGNAL')
	column = FPro_ColumnSeries(label, 1, ncolumn, 1, 1, 0, 0, 0,
     &			.TRUE.)
	star = FPro_LSplineFunc(column)
	call FPro_FreeObject(column)

c Set up linear splines for background signal
	ncolumn = FPro_LabelFind(label, 1, 'BACKGROUND_SIGNAL')
	column = FPro_ColumnSeries(label, 1, ncolumn, 1, 1, 0, 0, 0,
     &			.TRUE.)
	back = FPro_LSplineFunc(column)
	call FPro_FreeObject(column)

c Close file
	call FPro_FreeObject(label)

c Update domain
	xsize = FPro_ObjectOverlap(star, series, x1, x2)
	xsize = FPro_ObjectOverlap(back, series, x1, x2)

c***********************************************************************
c Get radius limits
c***********************************************************************

c Window curve
	curve2 = FPro_WindowFunc(curve, x1, x2)
	call FPro_FreeObject(curve)
	curve = curve2

c Summarize curve segments
	nsegments = FPro_CurveSegments(curve)
	do 200 nsegment = 0, nsegments
	    x = FPro_CurveExtremum(curve, nsegment, r, type)
	    y = FPro_CurveValue(curve, x, drdx)

	    write(*,20) nsegment, r, drdx*dx
20	    format(1x, i2, ': r = ', f10.3, ', dr = ', f8.3)
200	continue

c Get sampling parameters
	if (nsegments .gt. 1) then
	    write(*,10) 'Segment index = '
	    read(*,*,end=9999) nsegment
	else
	    nsegment = 1
	endif

	write(*,10) 'Radius min, max = '
	read(*,*,end=9999) r1, r2

	write(*,10) 'Resampling factor = '
	read(*,*,end=9999) expand

c***********************************************************************
c Set up filter
c***********************************************************************

c Set up new PSF
	write(*,10) 'Filter type (boxcar,triangle,sinc): (B/T/S) = '
	read(*,11,end=9999) char

	if (char .eq. 'B' .or. char .eq. 'b') then
	    newpsf = FPro_BoxcarFunc(1.d0, 0.5d0)
	else if (char .eq. 'T' .or. char .eq. 't') then
	    newpsf = FPro_TriangleFunc(1.d0, 1.d0)
	else if (char .eq. 'S' .or. char .eq. 's') then
	    newpsf = FPro_SincFunc(1.d0, 1.d0, SINC_HALFWIDTH, 0.d0)
	else
	    write(*,*) 'Unrecognized filter type: ' // char
	    goto 9999
	end if

c Old PSF for PPS data is a boxcar of 80% width
	oldpsf = FPro_BoxcarFunc(1.d0, 0.4d0)

c Center offset is 0. for odd expansion factors; 0.5 for even
	center = 0.d0
	if (expand .eq. 2*(expand/2)) center = 0.5d0

c Determine filter coefficients
	fsize = FPro_PSFFilter(oldpsf, newpsf, center, dble(expand),
     &		filter, NFILTER, f1, f2)

c We're done with PSF objects now
	call FPro_FreeObject(oldpsf)
	call FPro_FreeObject(newpsf)

c***********************************************************************
c Set up filtered data series
c***********************************************************************

c Now we're done with the original series
	call FPro_FreeObject(series)

c Create column object
c ntable=1, ncolumn=1, nitem=0 (all), row1=1, row2=0 (all), k1=0,
c buffersize=fsize, usedouble=FALSE
	column = FPro_ColumnSeries(dlabel, 1, 1, 0, 1, 0, 0, fsize,
     &			.FALSE.)

c Define uniformly weighted, uncorrelated samples
	series = FPro_FixedStat(column, 0, 1.d0)

c Set up filtered statistical series
	filtered = FPro_FilteredStat(series, filter, f1, f2)

c Free unneeded objects
	call FPro_FreeObject(dlabel)
	call FPro_FreeObject(column)
	call FPro_FreeObject(series)

c***********************************************************************
c Open output file
c***********************************************************************

	write(*,10) 'Output filename: '
	read(*,11,end=9999) filename

	open(OUTUNIT, file=filename, form='formatted', status='new'
     &		, carriagecontrol='list'
     &		)

c CARRIAGECONTROL is a useful qualifier supported by many
c implementations of FORTRAN.  Others seem to ignore it, but it may
c need to be commented out in some circumstances.

c***********************************************************************
c Convolve data at each selected point in array
c***********************************************************************

c Locate radius limits in sample coordinates
	xr1 = FPro_CurveInverse(curve, nsegment, r1)
	xr2 = FPro_CurveInverse(curve, nsegment, r2)

c Locate radius limits in series index
	kr1 = FPro_SeriesIndex(filtered, xr1) - center
	kr2 = FPro_SeriesIndex(filtered, xr2) - center

c Determine index direction and step
	dksign = 1
	if (kr2 .lt. kr1) dksign = -1

	dk = dksign * expand

c Shift index range inward to nearest integers
	k1 = nint(kr1)
	if ((k1 - kr1)*dksign .lt. 0.d0) k1 = k1 + dksign

	k2 = nint(kr2)
	if ((k2 - kr2)*dksign .gt. 0.d0) k2 = k2 - dksign

c Loop through filtered series, tabulating results
	do 300 k = k1, k2, dk
	    x = FPro_SeriesXValue(filtered, k+center)
	    r = FPro_FuncValue(curve, x)
	    b = FPro_FuncValue(back, x)
	    s = FPro_FuncValue(star, x)

	    count = FPro_SeriesValue(filtered, k, flag)
	    if (flag .eq. 0) then
		covar = FPro_StatCovar(filtered, k, k)
		sigma = sqrt(count * covar)

		tau    = TauFunc(count,         b, s, mu, MAXTAU)
		tau_lo = TauFunc(count + sigma, b, s, mu, MAXTAU)
		tau_hi = TauFunc(count - sigma, b, s, mu, MAXTAU)
	    else
		count  = 0.d0
		sigma  = MAXSIGMA

		tau    = 0.d0
		tau_lo = MINTAU
		tau_hi = MAXTAU
	    end if

	    write(OUTUNIT, 30) r, count, sigma,
     &	                       tau, tau_lo, tau_hi
30	    format(f9.2, ',', f7.3, ',', f6.3, ',',
     &			f7.4, ',', f7.4, ',', f7.4)
300	continue

	close(OUTUNIT)

c***********************************************************************
c Job is finished
c***********************************************************************

	call FPro_FreeObject(filtered)
	call FPro_FreeObject(curve)
	call FPro_FreeObject(back)
	call FPro_FreeObject(star)

9999	continue
	stop
	end

c***********************************************************************

c This function evaluates optical depth given the count value and
c calibration parameters.

	real*8 function TauFunc(count, b, s, mu, maxtau)
	real*8		count, b, s, mu, maxtau

	real*8		frac

	frac = (count - b) / s
	if (frac .le. 0.d0) then
	    TauFunc = maxtau
	    return
	end if

	TauFunc = -mu * log(frac)
	if (TauFunc .gt. maxtau) TauFunc = maxtau

	return
	end

c*******************************************************************************
c