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