c ppsresam.for
c*******************************************************************************
c$ Component_name:
c	PPSRESAM
c$ Abstract:
c	This program generates an ASCII-format file containing the PPS data
c	calibrated and resampled onto a uniform radius grid.
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 and interval (km).
c	(6) 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	(7) Output file name.
c$ Outputs:
c	A resampled 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 the PPS data
c	resampled onto a uniform radius grid.  The resultant file can be
c	plotted or used for further analysis.  This program also serves as a
c	sample program 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	RESAMOUT.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, interval (km) = 133000,134000,2
c	Filter type (boxcar,triangle,sinc): (B/T/S) = B
c	Output filename: RESAMOUT.TAB
c
c	The first five lines of the output file are as follows:
c	133000.00, 24.063, 1.052, 0.3891, 0.3608, 0.4191
c	133002.00, 24.601, 1.067, 0.3744, 0.3466, 0.4039
c	133004.00, 19.497, 0.958, 0.5357, 0.5010, 0.5731
c	133006.00, 21.232, 1.000, 0.4746, 0.4426, 0.5089
c	133008.00, 23.016, 1.032, 0.4189, 0.3894, 0.4503
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 Ring-Moon Systems 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 PPSRESAM
	implicit none
	include 'fprofile.inc'

	real*8		x, x1, x2, dx, xsize, y, r, r1, r2, dr, drdx,
     &			count, covar, sigma, b, s, mu, tau, tau_lo,
     &			tau_hi
	integer*4	ncolumn, nsegment, nsegments, type, fsize,
     &			nsamples, k, flag
	character	filename*80, char*1
	integer*4	dlabel, label, column, series, curve, curve2,
     &			star, back, newpsf, oldpsf, resampled
	real*8		TauFunc

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

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, interval (km) = '
	read(*,*,end=9999) r1, r2, dr

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***********************************************************************
c Set up resampled data series
c***********************************************************************

c Determine size of buffer
	fsize = FPro_ResampledSize(series, oldpsf, newpsf, curve,
     &			nsegment, r1, r2, dr)

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

c Create new column object with correct buffer size
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 resampled statistical series
	nsamples = (r2 - r1 + 1.01*dr) / dr
	resampled = FPro_ResampledStat(series, oldpsf, newpsf, curve,
     &			nsegment, r1, dr, 1, nsamples)

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

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 Write output file
c***********************************************************************

c Loop through resampled series
	do 300 k = 1, nsamples
	    r = FPro_SeriesXValue(resampled, dble(k))
	    x = FPro_CurveInverse(curve, nsegment, r)
	    b = FPro_FuncValue(back, x)
	    s = FPro_FuncValue(star, x)

	    count = FPro_SeriesValue(resampled, k, flag)

	    if (flag .eq. 0) then
		covar = FPro_StatCovar(resampled, 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(resampled)
	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