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 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 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