Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.atnf.csiro.au/computing/software/gipsy/tsk/velsmo.shl Дата изменения: Fri Jul 21 11:12:22 2000 Дата индексирования: Fri Jan 16 01:10:47 2009 Кодировка: Поисковые слова: aurora

E PROGRAM VELSMO (copyright notice)
C velsmo.shl
C COPYRIGHT (c) 1990
C Kapteyn Astronomical Institute - University of Groningen
C P.O. box 800, 9700 AV Groningen, The Netherlands
C
E PROGRAM VELSMO (velsmo.dc1)
C#> velsmo.dc1
C
CProgram: VELSMO
C
CPurpose: Program does a one dimensional smoothing of subsets.
C
CCategory: MANIPULATION
C
CFile: velsmo.shl
C
CAuthor: K.G. Begeman
C
CKeywords:
C
C INSET= Set and subsets to smooth. Maximum number of input and
C output subsets is 2048.
C
C BOX= Frame of output subsets [input subsets]. This will be
C the size of the output subsets.
C
C WEIGHTS= Give weights for neighbouring subsets [0.5].
C The first weight is for the closest neighbour subsets,
C the second weight for the second closest neighbour subsets
C etc. The weights will be normalized automatically, unless
C NORMALIZE=N is given. The maximum number of weights is
C the minimum of 100 and (number of subsets-1)/2.
C
C** NORMALIZE= Normalize the convolution function [Y]?
C
C OUTSET= Set and subsets where result is stored. Note that the
C number of output subsets is equal to the number of
C input subsets minus twice the number of weights entered.
C
CExample: VELSMO
C VELSMO Version 1.0 Feb 25, 1990
C VELSMO ,INSET=AURORA FREQ
C Set AURORA has 3 axes
C RA-NCP from -7 to 8
C DEC-NCP from -7 to 8
C FREQ-OHEL from 1 to 59
C VELSMO ,BOX=PC D 8 8
C BOX range for set AURORA :
C RA-NCP from -3 to 4
C DEC-NCP from -3 to 4
C VELSMO ,WEIGHTS=
C Convolution function: 0.25 0.50 0.25
C VELSMO ,OUTSET=HANNING
C Set HANNING has 3 axes
C RA-NCP from -3 to 4
C DEC-NCP from -3 to 4
C FREQ-OHEL from 2 to 58
C VELSMO +++ FINISHED +++
C
CUpdates: Feb 25, 1990: KGB, Document created.
C Dec 15, 1993: KGB, Keyword NOARMALIZE added.
C Jul 23, 1999: VOG, maxfie increased from 5 to 10
C Jul 14, 2000: VOG, maxfie increased from 10 to 256
C
C#<
E PROGRAM VELSMO (code)
program velsmo
C
C Declaration of parameters:
C
character*(*) ident
N Change version number on this line
parameter (ident = ' VELSMO Version 1.0 Feb 25, 1990 ')
integer maxaxes
N Maximum number of axes in a set
parameter (maxaxes = 10)
integer maxsubs
N Maximum number of subsets
parameter (maxsubs = 2048)
integer maxbuf
N Size of input/output data array
parameter (maxbuf = 4096)
integer maxfie
N Size of convolution function
parameter (maxfie = 256)
C
C Declarations for input set
C
N Name of input set
character*80 set1
N Permutation of axes for input set
integer axperm1(maxaxes)
N Ax size array
integer axsize1(maxaxes)
N Coordinate words for gdsi_read
integer cwlo1(maxsubs), cwup1(maxsubs)
N Frame of input subsets
integer blo(maxaxes), bup(maxaxes)
N Number of input subsets
integer nsub1
N Dimension of input set
integer setdim1
N Coordinate words of input subsets
integer subset1(maxsubs)
N Transfer id input set
integer tid1(maxsubs)
N Arrays for input data
real data1(maxbuf,-maxfie:maxfie)
C
C Declarations for output set
C
N Name of output set
character*80 set2
N Permutation of axes for output set
integer axperm2(maxaxes)
N Ax size array
integer axsize2(maxaxes)
N Coordinate words for gdsi_write
integer cwlo2(maxsubs), cwup2(maxsubs)
N Number of output subsets
integer nsub2
N Coordinate words of output subsets
integer subset2(maxsubs)
N Transfer id output set
integer tid2(maxsubs)
N Arrays for output data
real data2(maxbuf)
C
C Declaration of other variables:
C
N For formatted text output
character*512 text
N Circular addressing buffer
integer cirbuf(-maxfie:maxfie)
N Error return from various GDS routines
integer gerror
N Loop counter
integer m
N Counter for minmax3
integer mcount(maxsubs)
N Counter
integer n
N Number of blanks in output subset
integer nblank(maxsubs)
N Total number of pixels done of subset
integer ndone
N Length of convolution function
integer nfie
N Number of pixels read in
integer nread
N Counter for circular buffer
integer nrot
N Subset counters
integer ns1, ns2
N Total number of pixels in output subset
integer ntotal
N Dimension of subsets
integer subdim
N Normalization of convolution function
logical norm
N BLANK
real blank
N Convolution function
real confie(-maxfie:maxfie)
N Minimum and maximum in output subset
real datamin(maxsubs), datamax(maxsubs)
N Temps
real d1, d2
N Temp for convolution factor
real f
C
C Declaration of functions:
C
N Returns coordinate word
integer gdsc_fill
N Returns number of dimensions
integer gdsc_ndims
N Returns number of input subsets
integer gdsinp
N Returns number of output subsets
integer gdsout
N Returns number of logicals entered by user
integer userlog
N Returns number of reals entered by user
integer userreal
C
C Data statements:
C
N Error return from gds routines
data gerror / 0 /
N Set counters for minmax3
data mcount / maxsubs * 0.0 /
N Number of pixels done sofar
data ndone / 0 /
N Do normalize the convolution function
data norm / .TRUE. /
N Number of pixels in subset
data ntotal / 1 /
N Set transfer identifiers input data
data tid1 / maxsubs * 0 /
N Set transfer identifiers output data
data tid2 / maxsubs * 0 /
C
C Executable code:
C
N Get in touch with HERMES
call init
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Here we show the non-experienced user the version number and date
C of this program. Version 1.0 means that it has been migrated to
C portable GIPSY.
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N Tell user who we are
call anyout( 11, ident )
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Here we prompt the user for the name of the set (and subsets).
C The default message is used. VELSMO is a class 2 application so it
C will return the number of axis outside the subset.
C GDSINP will store the coordinate system of the chosen input set
C internally, so that GDSOUT can use it to create the output set.
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N Get number of input subsets
nsub1 = gdsinp( set1,
# subset1,
# maxsubs,
# 0,
# 'INSET=',
# ' ',
# 11,
# axperm1,
# axsize1,
# maxaxes,
# 2,
# 1 )
N Get dimension of input set
setdim1 = gdsc_ndims( set1, 0 )
N Get dimension of input subsets
subdim = gdsc_ndims( set1, subset1 )
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Get here the frame of the new subsets. Note that the size of
C the output subsets may be different.
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N Get subset frame
call gdsbox( blo,
# bup,
# set1,
# subset1,
# 1,
# 'BOX=',
# 'Frame of output subsets [input subset size]',
# 11,
# 0 )
N Find number of points in output subset
for n = 1, subdim
ntotal = ntotal * ( bup(n) - blo(n) + 1 )
cfor
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Now we calculate the maximum length of the convolution function.
C If only two input subsets were given, the programme will just
C transfer the subsets. The default convolution function is the
C HANNING taper (0.25, 0.50, 0.25 normalized).
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N Maximum length of convolution function
nfie = min( maxfie, ( nsub1 - 1 ) / 2 )
N Width of convolution function > one?
if (nfie .gt. 0)
then
N user request
nfie = userreal( confie(1),
# nfie,
# 1,
# 'WEIGHTS=',
# 'Give relative weights of neighbours [0.5]' )
N Default used?
if (nfie .eq. 0)
then
N Default weights = HANNING
confie(1) = 0.5
N Size of convolution function
nfie = 1
cif
cif
N Central weight
confie(0) = 1.0
N Number of output subsets
nsub2 = nsub1 - 2 * nfie
N Normalize?
n = userlog( norm,
# 1,
# 2,
# 'NORMALIZE=',
# 'Normalize convolution function [Y]?' )
if ( norm )
then
N Integrate the convolution function
for n = 1, nfie
confie(0) = confie(0) + 2.0 * confie(n)
cfor
N Normalize the convolution function
for n = 1, nfie
confie(n) = confie(n) / confie(0)
confie(-n) = confie(n)
cfor
N Normalize central weight
confie(0) = 1.0 / confie(0)
else
for n = 1, nfie
confie(-n) = confie(n)
cfor
cif
N Display convolution function
if (nfie .lt. 6)
then
write( text, '(''Convolution function:'',1000(f6.2:))')
# ( confie(n), n = -nfie, nfie )
N Send it out
call anyout( 11, text )
else
write( text, '(''Convolution function: ...'',
# 1000(f6.2:),''....'')')
# ( confie(n), n = -5, 5 )
N Send it out
call anyout( 11, text )
cif
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Here the internal buffer created by GDSINP is copied to a buffer
C which will be used by GDSOUT to create the output set (with
C exactly the same coordinate system as the input set, if the user
C wants to).
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N Assign GDSINP buffer to GDSOUT
call gdsasn( 'INSET=', 'OUTSET=', 1 )
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Now the subset size of the internal buffer for GDSOUT is changed
C according to what the user entered at the boxinp prompt. Note that
C only NAXIS% and CRPIX% need to be changed. The coordinate system
C is left completely unchanged.
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N Define size of output subsets
call gdscss( 'OUTSET=', blo, bup )
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Now we change the axis selection for the output set.
C ------------------------------------------------------------------
N Only last axis will be changed
axsize1(setdim1) = nsub2
N Change axis selection for output set
call gdscas( 'OUTSET=', subset1(nfie+1), axsize1 )
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C GDSOUT will create the output set if it does not exist, else it
C will only check if the subset size is exactly equal to the size
C specified with GDSCSS. In principal the user only has to give a
C set name, the rest is taken from the contents in the internal
C buffer of GDSOUT. If the user specifies also the subsets, GDSOUT
C will check whether they are present in the existing set or in the
C internal buffer if the subset does not exist.
C GDSOUT also copies the 'relevant' header information from the
C input set to the output set.
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N get output set
nsub2 = gdsout( set2,
# subset2,
# nsub2,
# 0,
# 'OUTSET=',
# ' ',
# 11,
# axperm2,
# axsize2,
# maxaxes )
N Loop which sets the input coordinate words
for ns1 = 1, nsub1
N Lower coordinate word of frame
cwlo1(ns1) = gdsc_fill( set1, subset1(ns1), blo )
N Upper coordinate word of frame
cwup1(ns1) = gdsc_fill( set1, subset1(ns1), bup )
cfor
N Loop which sets the output coordinate words
for ns2 = 1, nsub2
N Lower coordinate word of frame
cwlo2(ns2) = gdsc_fill( set2, subset2(ns2), blo )
N Upper coordinate word of frame
cwup2(ns2) = gdsc_fill( set2, subset2(ns2), bup )
cfor
N Import system defined BLANK
call setfblank( blank )
N Progress Bar
call stabar( 0.0, float(ntotal), float(ndone) )
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Here we repeat the operation for each subset
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
N loop
repeat
N reset input subset counter
ns1 = 1
N Read in enough data to produce one output subset
for n = -nfie, nfie
N Read data
call gdsi_read( set1,
# cwlo1(ns1),
# cwup1(ns1),
# data1(1,n),
# maxbuf,
# nread,
# tid1(ns1) )
N initialize circular buffer
cirbuf(n) = n
N Next input subset
ns1 = ns1 + 1
cfor
N Loop over output subsets
for ns2 = 1, nsub2
N Show Progress
call stabar( 0.0, float(ntotal),
# float(ndone) + float(nread*ns2) / float(nsub2) )
N Clear output buffer
call presetr( 0.0, data2, nread )
N Convolving loop
for n = -nfie, nfie
N Current factor
f = confie(n)
N Address of data array
l = cirbuf(n)
N Data loop
for m = 1, nread
N Temporary storage
d1 = data1(m,l)
d2 = data2(m)
N Already blank ?
if (d2 .ne. blank)
then
N Can we use it ?
if (d1 .ne. blank)
then
N Partial convolution
data2(m) = d2 + f * d1
else
N Blank forever
data2(m) = blank
cif
cif
cfor
cfor
N Determine minimum and maximum
call minmax3( data2,
# nread,
# datamin(ns2),
# datamax(ns2),
# nblank(ns2),
# mcount(ns2) )
N Write convolved data to disk
call gdsi_write( set2,
# cwlo2(ns2),
# cwup2(ns2),
# data2,
# nread,
# nread,
# tid2(ns2) )
N Rotate circular buffer
for n = -nfie, nfie
nrot = cirbuf(n) + 1
if (nrot .gt. nfie)
then
nrot = -nfie
cif
cirbuf(n) = nrot
cfor
N Read in data from next input subset
if (ns1 .le. nsub1)
then
N Read
call gdsi_read( set1,
# cwlo1(ns1),
# cwup1(ns1),
# data1(1,cirbuf(nfie)),
# maxbuf,
# nread,
# tid1(ns1) )
N Next subset
ns1 = ns1 + 1
cif
cfor
N New number of pixels done
ndone = ndone + nread
until (ndone .eq. ntotal)
N Update minimum, maximum and number of blank
call wminmax( set2,
# subset2,
# datamin,
# datamax,
# nblank,
# nsub2,
# 1 )
N Tell HERMES we're done
call finis
C
C End of program
C
stop
end