Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.atnf.csiro.au/computing/software/gipsy/tsk/combin.shl Дата изменения: Tue Feb 21 18:35:08 1995 Дата индексирования: Fri Jan 16 01:11:44 2009 Кодировка: Поисковые слова: aurora

C
C COPYRIGHT (c) 1990
C Kapteyn Astronomical Institute
C University of Groningen - 9700 AV Groningen, The Netherlands
C
C#> combin.dc1
C
CProgram: COMBIN
C
CPurpose: Combine sets to create new sets using mathematical
C expressions.
C
CCategory: ANALYSIS, COMBINATION, MANIPULATION, CALCULATION
C
CFile: combin.shl
C
CAuthor: M. Vogelaar
C
CKeywords:
C
C RESULT01= Combination code for first output set.
C The code must contain at least one parameter ($n).
C If you want a set with for example only noise, try
C something like: ($1 - $1) + expression for noise.
C
C RESULT02= Combination code for the second output set.
C If carriage return is given, COMBIN assumes that
C no other output combinations are wanted.
C
C RESULT03= etc.
C
C SET01= Input set (and subsets) for first combination.
C Maximum number of subsets is 2048. The output
C set(s) will get the same sizes as this first
C input set.
C
C BOX01= Frame for first input parameter. [entire subset]
C
C SET02= Input set (and subsets) for second parameter.
C Input is accepted if the number of subsets is
C equal to the previous number of subsets.
C If you want to use COMBIN outside a given box
C (see OPTION= keyword), all the input sets need
C to have the same axis sizes.
C
C BOX02= Frame for second parameter. [previous box]
C The axis lengths of this box have to be equal to the
C axis lengths of the previous box. If using COMBIN
C outside a box, the relative position of this box
C has to be the same as the position of the first box.
C
C SET03= etc.
C
C BOX03= etc.
C
C SETOUT01= Give set and subset(s) where the first combination is
C to be stored.
C
C SETOUT02= Give set and subset(s) where the second combination is
C to be stored.
C
C SETOUT03= etc.
C
C*** OPTION= [0]
C 0. COMBIN operates inside given box and blanks values
C outside this box
C 1. COMBIN operates inside given box and transfers
C 'outside' values of SET01
C 2. COMBIN operates outside given box and blanks values
C inside this box
C 3. COMBIN operates outside given box and transfers
C 'inside' values of SET01
C
C
CDescription: This program allows the user to create sets which are
C a (complicated) combination of other sets. The
C combination description language is best illustrated
C by the following example: suppose that one
C has two sets which are the result of a Fourier
C transform, one containing the real part and the other
C the imaginary part of the result. The code needed to
C convert these two sets to amplitude and phase is the
C following:
C
C RESULT01= SQRT( $1*$1 + $2*$2 ) (amplitude)
C RESULT02= DEG( ATAN2 ( $2, $1 ) ) (phase in degrees)
C
C Here $1 and $2 denote the first and second input sets.
C In this example two input sets result in two output
C sets. But just as easy we could produce a third output
C set which contains the scaled real part. Then the code
C would be:
C
C RESULT03= 2 * $1 + 1 (scaled real part)
C
C So with COMBIN we can combine M input sets into N
C output sets.
C The example above can be achieved with the following
C keywords:
C
C COMBIN,RESULT01=sqrt($1*$1+$2*$2)
C COMBIN,RESULT02=deg(atan2($2,$1))
C COMBIN,RESULT03=
C
C Note: Two different types of output sets will be
C produced
C
C Input sets:
C
C COMBIN,SET01=AURORA_RE FREQ 1:10
C COMBIN,SET02=AURORA_IM FREQ 1:10
C
C Output sets:
C
C COMBIN,SETOUT01=AURORA_AM FREQ 1:10
C COMBIN,SETOUT02=AURORA_PH FREQ 1:10
C
C Note: We used (the default) option 0
C
C Available mathematics:
C
C 1) functions; syntax ff(..); where ff is one of the
C following available functions:
C abs(x) absolute value of x
C sqrt(x) square root of x
C sin(x) sine of x
C asin(x) inverse sine of x
C cos(x) cosine of x
C acos(x) inverse cosine of x
C tan(x) tangent of x
C atan(x) inverse tan of x
C atan2(x,y) inverse tan (mod 2pi) x = sin, y = cos
C exp(x) exponential of x
C ln(x) natural log of x
C log(x) log (base 10) of x
C sinh(x) hyperbolic sine of x
C cosh(x) hyperbolic cosine of x
C tanh(x) hyperbolic tangent of x
C rad(x) convert x to radians
C deg(x) convert x to degrees
C erf(x) error function of x
C erfc(x) 1-error function
C max(x,y) maximum of x and y
C min(x,y) minimum of x and y
C sinc(x) sin(x)/x (sinc function)
C sign(x) sign of x (-1,0,1)
C mod(x,y) gives remainder of x/y
C int(x) truncates to integer
C nint(x) nearest integer
C ranu(x,y) generates uniform noise between x and y
C rang(x,y) generates gaussian noise with mean x
C and dispersion y
C ranp(x) generates poisson noise with mean x
C ifeq(x,y,a,b) returns a if x == y, else b
C ifne(x,y,a,b) returns a if x != y, else b
C ifgt(x,y,a,b) returns a if x > y, else b
C ifge(x,y,a,b) returns a if x >= y, else b
C iflt(x,y,a,b) returns a if x < y, else b
C ifle(x,y,a,b) returns a if x <= y, else b
C ifblank(x,a,b) returns a if x == BLANK, else b
C
C Some (statistical) functions have a variable number
C of arguments:
C
C sum(x1,..,xn) sum of elements x1 .. xn
C mean(x1,..,xn) mean
C var(x1,..,xn) variance
C sdev(x1,..,xn) standard deviation
C adev(x1,..,xn) absolute deviation
C skew(x1,..,xn) skewness
C kurt(x1,..,xn) kurtosis
C median(x1,..,xn) median
C nblanks(x1,..,xn) number of blanks
C
C max(x1,..,xn) maximum of elements x1 .. xn
C min(x1,..,xn) minimum
C
C Note that n <= 128
C
C 2) constants; syntax cc; where cc is one of the following
C available constants:
C PI 3.14159....
C C speed of light (SI)
C H Planck (SI)
C K Boltzmann (SI)
C G gravitation (SI)
C S Stefan-Boltzman (SI)
C M mass of sun (SI)
C P parsec (SI)
C BLANK Universal undefined value
C Note: the Hubble constant is not included.
C 3) operators; syntax op; where op is one of the following
C available operators:
C + addition
C - subtraction
C * multiplication
C / division
C ** power
C 4) parameters; syntax $n; where 1 <= n <= 128.
C
CRemarks: A) The calculations are all done in double precision.
C
C B) BLANK values are recognized. Any operation on a
C BLANK causes the result to be set to BLANK (except the
C function IFBLANK). Also all illegal operations, like
C dividing by zero or the logarithm of a negative value,
C will cause the result to be set to BLANK.
C
C C) If you cannot find your favorite constant or function
C in the list, please contact Kor Begeman. He might be
C persuaded to put it in.
C
CExamples: 1) Add three maps
C RESULT01=$1+$2+$3
C RESULT02=
C SET01= set and subset(s) of first parameter
C SET02= set and subset(s) of second parameter
C SET03= set and subset(s) of third parameter
C SETOUT01= set and subset(s) of sum
C
C 2) Add gaussian noise to a map with zero mean and a
C rms of 2
C RESULT01=$1+rang(0,2)
C RESULT02=
C SET01= input set and subset(s)
C SETOUT01= output set and subset(s) with noise added
C
C 3) Clip values in a set between -4 and 5 and replace with
C 'BLANK'
C RESULT01=ifgt($1,5,BLANK,iflt($1,-4,BLANK,$1))
C RESULT02=
C SET01= input set and subset(s)
C SETOUT01= clipped output set and subset(s) with values
C between -4 and 5
C
C 4) Conditional transfer: if value in second set is
C greater than 2 then result is value from first set,
C else set this value to BLANK
C RESULT01=ifgt($2,2,$1,BLANK)
C RESULT02=
C SET01= input set and subset(s)
C SET02= test set and subset(s)
C SETOUT01= conditional transfered output set and
C subset(s)
C
C 5) Calculate median of three sets and put result in output
C set MEDIAN
C RESULT01=median($1,$2,$3)
C RESULT02=
C SET01=OPTC1
C SET01=OPTC2
C SET01=OPTC3
C SETOUT01=MEDIAN
C
CNotes:
C#<
C
E ***** Program COMBIN *****
PROGRAM COMBIN

C Declaration of parameters:

CHARACTER*(*) ident
PARAMETER ( ident = ' Version 1.0 Mar 14, 1990 ' )
INTEGER maxsubs
PARAMETER ( maxsubs = 2048 )
INTEGER maxaxes
PARAMETER ( maxaxes = 10 )
N Default options for use in USERxxx() routines
INTEGER none, request, hidden
PARAMETER ( none = 0 )
PARAMETER ( request = 1 )
PARAMETER ( hidden = 2 )
N Buffer size for I/O
INTEGER maxIObuf
PARAMETER ( maxIObuf = 4096 )
N Remove with 'WMINMAX' old minmax descriptors
INTEGER remove
PARAMETER ( remove = 1 )
N Max. number of characters in 1 comb. description
INTEGER maxchars
PARAMETER ( maxchars = 256 )
N 'FIE' has 16 different buffers for results
INTEGER maxrslt
PARAMETER ( maxrslt = 16 )
N Syntax parameters: $n with 1 <= n <= 32
INTEGER maxpars
PARAMETER ( maxpars = 32 )
N Coordinate word for whole set
INTEGER wholeset
PARAMETER ( wholeset = 0 )

C Declarations for GDSINP:

INTEGER GDSINP
CHARACTER*80 setname( maxpars )
N Array containing subset coordinate words
INTEGER subsets( maxsubs, maxpars )
N Axis permutation array
INTEGER axperm( maxaxes, maxpars )
INTEGER axcount( maxaxes, maxpars )
N Number of subsets
INTEGER nsubs
N Determined at first input of sets
INTEGER nsubsmax
INTEGER dfault
N Number of output device [0..16]
INTEGER devicenum
N Keywords for numerous inputs
CHARACTER*10 keyword
CHARACTER*10 setkwrd
CHARACTER*10 boxkwrd
CHARACTER*40 message
N Algorithm with repeat axes ==> class = 1
INTEGER class
N Dimension of the subsets
INTEGER subdim

C Declarations for GDSBOX:

N Grid vectors of sub frames of input sets
INTEGER BgridLO( maxaxes, maxpars )
INTEGER BgridHI( maxaxes, maxpars )
N Grid vectors of entire frames of inp. sets
INTEGER FgridLO( maxaxes, maxpars )
INTEGER FgridHI( maxaxes, maxpars )
N Same for output sets
INTEGER FgridLOO( maxaxes, maxrslt )
INTEGER FgridHIO( maxaxes, maxrslt )
N What's the default in GDSBOX?
INTEGER boxopt

C Declarations for GDSOUT:

INTEGER GDSOUT
N Number of output sets
INTEGER nsubsO
INTEGER subsetsO( maxsubs, maxrslt )
INTEGER axpermO( maxaxes, maxrslt )
INTEGER axcountO( maxaxes, maxrslt )
CHARACTER*80 setnameO( maxrslt )

C Functions:

N Extracts grid coordinate from coord. word
INTEGER GDSC_GRID
N Returns coordinate word
INTEGER GDSC_FILL
N Returns .TRUE. if inside defined sub frame
LOGICAL OUTSIDEPTR
N Returns .TRUE. if outside defined sub frame
LOGICAL INSIDEPTR
N Determine length of string
INTEGER NELC
N Obtain the name under which a GIPSY task runs
CHARACTER*9 MYNAME
N Get length of input string
INTEGER USERTEXT
N Decode a string with math expression for FIEDO
INTEGER FIEINI
N Evaluate the code generated by FIEINI
INTEGER FIEDO
N User integer input
INTEGER USERINT

C Variables for the algorithms:

N Various counters
INTEGER m, i, q
N There are totpixels in one subset
INTEGER totpixels
N Array version of transfer id's
INTEGER tids( maxpars )
N For output data
INTEGER tidsO( maxrslt )
N Special tid for outside operations
INTEGER tidT
N Counter for subroutine MINMAX3
INTEGER mcount( maxrslt )
N Number of pointers in INITPTR buffer
INTEGER ptrcount
INTEGER stilltowrite
N Coordinate words for total & sub frame
INTEGER Fcwlo( maxpars )
INTEGER Fcwhi( maxpars )
N Array version of Bcwlo, Bcwhi
INTEGER Bcwlo( maxpars )
INTEGER Bcwhi( maxpars )
N Array version of FcwloO, FcwhiO
INTEGER FcwloO( maxrslt )
INTEGER FcwhiO( maxrslt )
N For use in INSIDEPTR/OUTSIDEPTR
INTEGER bufptr
N Dynamical I/O buffer length
INTEGER numinreadbuf
N Number of elements inside/outside sub frame
INTEGER numpixin
INTEGER numpixout
N Number of pixels actually in read buffer
INTEGER numpixread
N How far we proceeded, used in STABAR
REAL curval
N Help variable in STABAR
REAL ratio
N Offset for 1 dimensional dataI array
INTEGER offset
N For updating number of blanks
INTEGER nblanks( maxsubs, maxrslt )
N For updating minimum and maximum of subset
REAL minval( maxsubs, maxrslt )
REAL maxval( maxsubs, maxrslt )
N Arrays holding the data
REAL dataO( maxIObuf, maxrslt )
REAL dataI( maxIObuf * maxpars )
N Data array only for transfer option
REAL dataT( maxIObuf )
N Status of function FIEDO
INTEGER fieresult

C Miscellaneous:

N Error codes for Range and Grid routines
INTEGER Rerror, Gerror
N Coordinate words to determine total range
INTEGER cwlo, cwhi
N Option 0..3, for different operations
INTEGER option
N Choose COMBIN operating in- or outside box
LOGICAL workinside
N Choose between transfer or blanking
LOGICAL transfer
N True when grids of input edges are equal
LOGICAL equalgrids
N Loop guard for valid set input
LOGICAL agreed
N Loop guard for valid box input
LOGICAL equalbox
N Loop guard for 'work outside' option
LOGICAL equalframe
N Another guard for 'work outside' option
LOGICAL equalpos
N First input of a valid set
LOGICAL first
N Length of result string
INTEGER len
N Number of results (=output sets)
INTEGER numouts
N Counter for 'numouts'
INTEGER out
N Total number of output sets
INTEGER totouts
N Number of parameters = num. of input sets
INTEGER numpars
N Maximum number of counted parameters
INTEGER totpars
N Id. for each result string
INTEGER fid( maxrslt )
N Position of error in result string
INTEGER errpos
N Current parameter
INTEGER par
N Length of axes of first box
INTEGER Blen1( maxpars )
N Length of one axis of other boxes
INTEGER Blen2
N Display name of this program
CHARACTER*9 task
N Storage for 'maxrslt' descriptions
CHARACTER*( maxchars ) resultstr( maxrslt )


E ***** COMBIN: Main *****

N Enter HERMES
CALL INIT
N Get task name
task = MYNAME()
N Show task and version number
CALL ANYOUT( 8, ( task(:NELC(task) ) // ident) )

C Options for work area, asked before(!) GDSBOX:

keyword ='OPTION='
message ='Give option: [0]'
dfault = hidden

C Choose one of the four available options:

REPEAT
IF ( USERINT( option, 1, dfault, keyword, message ) .EQ. 0 )
THEN
option = 0
CIF
agreed = ( (option .GE. 0) .AND. (option .LE. 3) )
IF (.NOT. agreed)
THEN
dfault = request
message = 'COMBIN: Illegal opt. try again: [0]'
CALL CANCEL( keyword )
CIF
UNTIL agreed

C Translate options into logicals to work with:

IF (option .EQ. 0)
THEN
workinside = .TRUE.
transfer = .FALSE.
CALL ANYOUT( 3, 'Option: operate inside box, blank outside' )
CIF
IF (option .EQ. 1)
THEN
workinside = .TRUE.
transfer = .TRUE.
CALL ANYOUT( 3, 'Option: operate inside box, transfer outside' )
CIF
IF (option .EQ. 2)
THEN
workinside = .FALSE.
transfer = .FALSE.
CALL ANYOUT(3, 'Option: operate outside box, blank inside' )
CALL ANYOUT(3, 'Note: axis sizes of frames need to be equal!')
CIF
IF (option .EQ. 3)
THEN
workinside = .FALSE.
transfer = .TRUE.
CALL ANYOUT(3, 'Option: operate outside box, transfer inside' )
CALL ANYOUT(3, 'Note: axis sizes of frames need to be equal!')
CIF

C Input of combination descriptions:

numouts = 1
totpars = 0
keyword = 'RESULT01='
message = 'Give expression: '
dfault = none
REPEAT
fid(numouts) = 0
resultstr(numouts) = ' '
len = USERTEXT( resultstr(numouts), dfault, keyword, message )
N Check string & prepare for next input
IF (len .NE. 0)
THEN
N Check whether input string is ok
numpars = FIEINI( resultstr(numouts), fid(numouts),
# errpos )
N Store the maximum value of numpars in totpars
totpars = MAX( numpars, totpars )
IF (numpars .GT. 0)
THEN
numouts = numouts + 1
dfault = request
ELSE
CALL CANCEL( keyword )
IF (numpars .EQ. 0)
THEN
CALL ANYOUT( 3, 'COMBIN: No parameters specified' )
CIF
IF (numpars .LT. 0)
THEN
WRITE( message,
# '(''COMBIN: Error in string at position '', I3 )' )
# errpos
CALL ANYOUT( 3, message )
CIF
CIF
WRITE( keyword, '(''RESULT'', I2.2, ''='')' ) numouts
WRITE( message, '(''Give expression: '', I2, '' [stop]'' )' )
# numouts
CIF
UNTIL ( ( (numouts-1) .EQ. maxrslt) .OR. (len .EQ. 0) )
N Compensate for previous loop
totouts = numouts - 1

C Get sets corresponding to the parameters:

N Prepare GDSINP for first SET
dfault = none
N Application repeats operation for each subset
class = 1
N First time free choice, later not!
subdim = 0
devicenum = 11
N Do it for all input parameters
FOR par = 1, totpars
first = (par .EQ. 1)
WRITE( setkwrd, '(''SET'', I2.2, ''='' )' ) par
WRITE( message, '(''Set (and subsets) of parameter '', I2 )' )
# par
N Until nsubs eq. to first time number of subs.
REPEAT
nsubs = GDSINP( setname(par), subsets(1, par), maxsubs,
# dfault, setkwrd, message, devicenum,
# axperm(1, par), axcount(1, par),
# maxaxes, class, subdim )
IF first
N Determine axis lengths & first # of subsets
THEN
nsubsmax = nsubs
totpixels = 1
FOR m = 1, subdim
totpixels = totpixels * axcount(m, 1)
CFOR
CIF
agreed = (nsubs .EQ. nsubsmax)
IF (.NOT. agreed)
THEN
CALL ANYOUT( 3, 'COMBIN: Wrong number of subsets' )
CALL CANCEL( setkwrd )
N Else if agreed get grids
ELSE
Rerror = 0
Gerror = 0
CALL GDSC_RANGE( setname(par), 0, cwlo, cwhi, Rerror )
FOR m = 1, subdim
FgridLO(m, par) = GDSC_GRID( setname(par),
# axperm(m, par), cwlo, Gerror )
FgridHI(m, par) = GDSC_GRID( setname(par),
# axperm(m, par), cwhi, Gerror )
CFOR
equalgrids = .TRUE.
IF (.NOT. first)
THEN
N Are frst and scnd sets comparable on sset level?
FOR q = 1, subdim
equalgrids = ( equalgrids .AND.
# ( FgridLO(q, par-1) .EQ.
# FgridLO(q, par)
# ) .AND.
# ( FgridHI(q, par-1) .EQ.
# FgridHI(q, par)
# ) )

CFOR
N For outside operations eq. Fsizes are needed
IF ( (.NOT. workinside) .AND. (.NOT. equalgrids) )
THEN
N Compare sizes of frames
equalframe = .TRUE.
FOR q = 1, subdim
equalframe = (equalframe .AND.
# ( axcount(q, par) .EQ. axcount(q, 1) ) )
CFOR
IF (.NOT. equalframe)
THEN
agreed = .FALSE.
CALL ANYOUT( 3, 'COMBIN: For this option you need' //
# ' same axis lengths as ' )
CALL ANYOUT( 3, ' previous set frame!' )
CALL CANCEL( setkwrd )
CIF
N End if work outside option
CIF
N End of not first
CIF
N End of agreed / not agreed
CIF
UNTIL agreed

C Input of sub frames with GDSBOX:

N Prepare for a frame for first set
dfault = request
WRITE( boxkwrd, '(''BOX'', I2.2, ''='' )' ) par
message = ' '
N Default is entire subset
IF first
THEN
boxopt = 0
CIF
devicenum = 11
REPEAT
N Set default sizes for GDSBOX
IF (.NOT. first)
THEN
N Fill arrays with GRIDS for GDSBOX
FOR q = 1, subdim
BgridHI(q, par) = BgridHI(q, par - 1)
BgridLO(q, par) = BgridLO(q, par - 1)
CFOR
N Default in B..LO & B..HI
boxopt = 6
IF (.NOT. equalgrids)
N Fill arrays with SIZES for GDSBOX
THEN
FOR q = 1, subdim
BgridHI(q,par) = BgridHI(q,par-1)-BgridLO(q,par-1) + 1
CFOR
N Box restricted to size in B..HI
boxopt = ( 8 + 4 )
N End of not equalgrids
CIF
N End if not first
CIF
CALL GDSBOX( BgridLO(1, par), BgridHI(1, par),
# setname(par), subsets(1, par),
# dfault, boxkwrd, message, devicenum, boxopt )
equalpos = .TRUE.
IF first
N Determine length of box axes of first box
THEN
equalbox = .TRUE.
FOR q = 1, subdim
Blen1(q) = BgridHI(q, 1) - BgridLO(q, 1) + 1
CFOR
ELSE
N Not the first ==> compare box sizes
equalbox = .TRUE.
FOR q = 1, subdim
Blen2 = BgridHI(q, par) - BgridLO(q, par) + 1
equalbox = ( equalbox .AND. (Blen1(q) .EQ. Blen2) )
CFOR
N For outside operations, box must have same
equalpos = .TRUE.
N relative position as previous box
IF (equalbox .AND. (.NOT. workinside) )
THEN
FOR q = 1, subdim
equalpos = ( equalpos .AND. (
# (BgridLO(q, par) - FgridLO(q, par)) .EQ.
# (BgridLO(q, 1 ) - FgridLO(q, 1 )) ) )
CFOR
IF (.NOT. equalpos)
THEN
CALL ANYOUT( 3, 'COMBIN: Wrong position of box' )
CALL CANCEL( boxkwrd )
CIF
CIF
N End first or other set
CIF
IF (.NOT. equalbox)
THEN
CALL ANYOUT( 3, 'COMBIN: Wrong size of box!' )
CALL CANCEL( boxkwrd )
CIF
equalbox = (equalbox .AND. equalpos)
N Now its ok with the boxes also
UNTIL equalbox
N Close after 'totpars'
CFOR

C Create output sets:

FOR out = 1, totouts
WRITE( keyword, '(''SETOUT'', I2.2, ''='' )' ) out
N Assign GDSINP buffer of SET01 to GDSOUT
CALL GDSASN( 'SET01=', keyword, class )
dfault = none
message ='Give output set (and subsets):'
devicenum = 11
N Until the correct number of SETOUT subsets
REPEAT
nsubsO = GDSOUT( setnameO(out), subsetsO(1, out),
# nsubsmax, dfault, keyword, message,
# devicenum, axpermO(1, out),
# axcountO(1, out), maxaxes )
agreed = (nsubsO .EQ. nsubsmax)
IF (.NOT. agreed)
THEN
CALL ANYOUT( 3, 'COMBIN: Wrong number of subsets' )
CALL CANCEL( keyword )
CIF
UNTIL agreed
Rerror = 0
Gerror = 0
CALL GDSC_RANGE( setnameO(out), wholeset, cwlo, cwhi, Rerror )
FOR q = 1, subdim
FgridLOO(q, out) = GDSC_GRID( setnameO, axpermO(q, out),
# cwlo, Gerror )
FgridHIO(q, out) = GDSC_GRID( setnameO, axpermO(q, out),
# cwhi, Gerror )
CFOR
N End for, for all setouts
CFOR

C Fill buffers and calculate:

N cw=coordinate word, lo=lower, hi=upper
N F=Frame, B=Box
FOR m = 1, nsubsmax
N Determine coordinate words for input sets
FOR par = 1, totpars
Fcwlo(par) = GDSC_FILL( setname(par), subsets(m, par),
# FgridLO(1, par) )
Fcwhi(par) = GDSC_FILL( setname(par), subsets(m, par),
# FgridHI(1, par) )
Bcwlo(par) = GDSC_FILL( setname(par), subsets(m, par),
# BgridLO(1, par) )
Bcwhi(par) = GDSC_FILL( setname(par), subsets(m, par),
# BgridHI(1, par) )
tids(par) = 0
CFOR
N Determine coordinate words for output sets
FOR out = 1, totouts
FcwloO(out) = GDSC_FILL( setnameO(out), subsetsO(m, out),
# FgridLOO(1, out) )
FcwhiO(out) = GDSC_FILL( setnameO(out), subsetsO(m, out),
# FgridHIO(1, out) )
N Reset transfer id for output data
tidsO(out) = 0
N Reset counter for subroutine MINMAX3
mcount(out) = 0
CFOR
ptrcount = 0
numpixread = 0
tidT = 0
stilltowrite = totpixels
REPEAT
N # per iteration cannot exceed maxIObuf
numinreadbuf = MIN( maxIObuf, stilltowrite )
N On exit ptrcount contains # ptrs in buf.
CALL INITPTR( FgridLO(1, 1), FgridHI(1, 1),
# BgridLO(1, 1), BgridHI(1, 1),
# subdim, numinreadbuf, ptrcount )
N Read full buff., all first set data is needed
IF transfer
THEN
CALL GDSI_READ( setname(1), Fcwlo(1), Fcwhi(1),
# dataT, numinreadbuf, numpixread, tidT )
CIF
IF workinside
THEN
WHILE ( OUTSIDEPTR( bufptr, numpixout) )
N Then set to blank
IF (.NOT. transfer)
THEN
FOR out = 1, totouts
CALL SETNFBLANK( dataO(bufptr+1, out), numpixout )
CFOR
N Else if transfer is needed ...
ELSE
FOR out = 1, totouts
CALL MOVER( dataT(bufptr+1), dataO(bufptr+1, out),
# numpixout )
CFOR
CIF
CWHILE
numpixread = 0
WHILE ( INSIDEPTR( bufptr, numpixin ) )
N Read data of all parameters in a 1 dim. array
offset = 1
FOR par = 1, totpars
CALL GDSI_READ( setname(par),
# Bcwlo(par), Bcwhi(par),
# dataI(offset), numpixin, numpixread,
# tids(par) )
offset = offset + numpixread
CFOR
FOR out = 1, totouts
N Do the calculations
fieresult = FIEDO( dataI, numpixin,
# dataO(bufptr+1, out), fid(out) )
IF ( fieresult .LT. 0 )
THEN
CALL ERROR( 4, 'COMBIN: No evaluation possible' )
CIF
CFOR
CWHILE
N If NOT workinside ==> operate outside box
ELSE
offset = 1
FOR par = 1, totpars
N Combin on all(!) data, blank later
CALL GDSI_READ( setname(par),
# Fcwlo(par), Fcwhi(par),
# dataI(offset), numinreadbuf, numpixread,
# tids(par) )
offset = offset + numpixread
CFOR
FOR out = 1, totouts
fieresult = FIEDO( dataI, numpixread,
# dataO(1, out), fid(out) )
CFOR
N Blank the inside data
WHILE ( INSIDEPTR( bufptr, numpixin ) )
N Then set to blank
IF (.NOT. transfer)
THEN
FOR out = 1, totouts
CALL SETNFBLANK( dataO(bufptr+1, out), numpixin )
CFOR
N Else if transfer is needed ...
ELSE
FOR out = 1, totouts
CALL MOVER( dataT(bufptr+1), dataO(bufptr+1, out),
# numpixin )
CFOR
CIF
CWHILE
CIF
FOR out = 1, totouts
N Write output data
CALL GDSI_WRITE( setnameO(out), FcwloO(out), FcwhiO(out),
# dataO(1, out), numinreadbuf, numinreadbuf,
# tidsO(out) )
N Find running min, max etc for this output
CALL MINMAX3( dataO(1, out), numinreadbuf, minval(m, out),
N set and subset and store it in arrays
# maxval(m, out), nblanks(m, out), mcount(out) )
CFOR
stilltowrite = stilltowrite - numinreadbuf
ratio = FLOAT( totpixels-stilltowrite ) / FLOAT(totpixels)
curval = ( FLOAT( m-1 ) + ratio ) / FLOAT( nsubsmax )
CALL STABAR( 0.0, 1.0, curval )
UNTIL (stilltowrite .EQ. 0)
CFOR
N min max probably changed ==> update min, max
FOR out = 1, totouts
CALL WMINMAX( setnameO(out), subsetsO(1, out),
# minval(1, out), maxval(1, out), nblanks(1, out),
# nsubsmax, remove )
CFOR

N Exit HERMES
CALL FINIS

STOP
END


E ***** COMBIN Technicalities *****

C Options in GDSBOX:
C
C 1 box may exceed subset size
C 2 default is in B..LO
C 4 default is in B..HI
C 8 box restricted to size defined in B..HI
C These codes work additive.

C The size of the input set will be copied to the output set.
C It is possible to work on a smaller area. Therefore we need
C the grid coordinates of the frame ('F') and
C the coordinates of the user defined box ('B').

C All output sets will get same coordinate system as first input set

C The order of the indices of the two dim. arrays is important,
C that is, if they are used in a subroutine call. The most rapidly
C changing index is always the first.