subroutine bdmset (wav,rflib,cl1,cl2,cr1,cr2,
rflibc,minch,maxch,ifeattype,ier)
implicit integer*4 (i-n)
#ccc name: bdmset
#ccc version date: January 29, 1990
#ccc author(s): Roger N. Clark
#ccc language: Ratfor
#ccc
#ccc short description:
#ccc This program computes non-changing library reference data
#ccc for the bandmp subroutine.
#ccc
#ccc algorithm description:
#ccc given continuum removed library spectrum,
#ccc add a constant to the spectrum and renormalize
#ccc to modify band depth such that there is a least
#ccc squares fit to the observed spectrum.
#ccc
#ccc THIS subroutine computes the non-changing
#ccc components of the least squares solution.
#ccc system requirements: none
#ccc subroutines called: none
#ccc
#ccc
#ccc parameter description:
#ccc INPUT:
#ccc rflib = reference library spectrum R*4 (cr2 elements)
#ccc cl1 = continuum point begin on left side of band I*4
#ccc cl2 = continuum point end on left side of band I*4
#ccc note: cl2 >= cl1 (checked)
#ccc cr1 = continuum point begin on right side of band I*4
#ccc note: cr1 > cl2 + 1 (checked)
#ccc cr2 = continuum point end on right side of band I*4
#ccc note: cr2 >= cr1 (checked)
#ccc note: cr2 also determines the max array sizes
#ccc OUTPUT:
#ccc rflibc = reference library spectrum, continuum
#ccc removed R*4 (cr2 elements)
#ccc minch = minimum channel in the reference library spectrum.
#ccc This is defined to be the band minimum.
#ccc maxch = maximum channel in the reference library spectrum.
#ccc This is defined to be the band maximum.
#ccc ifeattype = -1 is an emission feature
#ccc = 1 is an absorption band
#ccc ier = error detected in input:
#ccc = 0 no error
#ccc = 1 error
#ccc
#ccc common description: none
#ccc message files referenced: none
#ccc internal variables:
#ccc file description: none
#ccc user command lines: none
#ccc update information:
#ccc NOTES:
#ccc
integer*4 ttyout
integer*4 cl1,cl2,cr1,cr2,ier,minch,maxch,ifeattype
real*4 wav(cr2), rflib(cr2)
real*4 rflibc(cr2)
ttyout = 6
minch=0
maxch=0
#
# check input parameters are correct
#
ier = 0
if (cl1 > cl2) {
write (ttyout, 100)
100 format (' ERROR: left continuum segment point 1 is',
' greater than point 2')
ier = 1
}
if (cr1 > cr2) {
write (ttyout, 101)
101 format (' ERROR: right continuum segment point 1 is',
' greater than point 2')
ier = 1
}
if (cl2 + 1 > cr1 - 1) {
write (ttyout, 102)
102 format (' ERROR: left continuum segment point 2 + 1 ',
'channel is greater than ',/,
' right continuum segment point 1 - 1 channel',/,
' This means there are no channels for',
' the actual absorption band' )
ier = 1
}
if (ier == 1) return
#
# deleted point value:
#
delpt = -1.23e34
#
# compute continuum correction to ref spectrum
#
# avlr = average left side of continuum
# avrc = average right side of continuum
# avwlc = average wavelength left side continuum
# avwrc = average wavelength right side continuum
# first left continuum average:
avlc = 0.0
avwlc = 0.0
n = 0
do i = cl1, cl2 {
if (wav(i) == delpt || rflib(i) == delpt) next
avlc = avlc + rflib(i)
avwlc = avwlc + wav(i)
n = n +1
}
if (n < 1) {
write (ttyout,120)
120 format (' ERROR: all points in reference spectrum,',
' left continuum deleted.')
ier = 1
return
}
avlc = avlc / float(n)
avwlc = avwlc / float(n)
# now right continuum average:
avrc = 0.0
avwrc = 0.0
n = 0
do i = cr1, cr2 {
if (wav(i) == delpt || rflib(i) == delpt) next
avrc = avrc + rflib(i)
avwrc = avwrc + wav(i)
n = n +1
}
if (n < 1) {
write (ttyout,121)
121 format (' ERROR: all points in reference spectrum,',
' right continuum deleted.')
ier = 1
return
}
avrc = avrc / float(n)
avwrc = avwrc / float(n)
# now we have two pairs of x,y points to derive the continuum line
# thus: wav = a * rflib + b
bottom = avwrc - avwlc
if (abs(bottom) < 0.1e-20) {
write (ttyout, 130)
130 format (' ERROR: wavelength range of continuum is',
' too small!')
ier = 1
return
}
a = (avrc - avlc)/bottom
b = avrc - a * avwrc
#
# now we can continuum correct the reference spectrum.
# compute rflibc
#
if (cl1 > 1) {
do i = 1, cl1-1 {
rflibc(i) = delpt
}
}
do i = cl1, cr2 {
if (wav(i) == delpt || rflib(i) == delpt) {
rflibc(i) = delpt
} else {
contin = a * wav(i) + b
if (abs(contin) > 0.1e-20) {
rflibc(i) = rflib(i) / contin
} else {
rflibc(i) = delpt
}
}
}
#
# now find the band minimum and maximum in case of emission feature
#
il = cl2 + 1
ir = cr1 - 1
195 if (rflibc(il) == -1.23e+34) { # find first non-deleted point
il = il + 1
if (il > ir) {
write (ttyout, 196)
196 format (' ERROR in band map setup:',/,
' can not find a min or max: ',
'all channels deleted')
ier = 1
return
}
go to 195
}
rmin = rflibc(il)
rmax = rmin
minch = il
maxch = il
do i = il, ir {
if (rflibc(i) == delpt) next
if (rflibc(i) < rmin) {
minch = i
rmin = rflibc(i)
}
if (rflibc(i) > rmax) {
maxch = i
rmax = rflibc(i)
}
}
depth = 1.0 - rmin # band depth for an emission feature
emiss = rmax - 1.0 # emission feature
if (rmin < 0.0) {
write (ttyout, 200)
200 format (' ERROR: band depth in the library spectrum ',/,
' is LESS THAN ZERO: that is invalid')
ier = 1
return
}
#
# determine if the feature is emission or absorption
#
ifeattype = 0
if (emiss > depth) {
write (ttyout, 201)
201 format (' NOTE: this feature is an emission feature.')
ifeattype = -1 # emission feature
} else {
ifeattype = 1 # absorption band
}
if (ifeattype == -1 & maxch == 0) {
write (ttyout,303) maxch
303 format ('ERROR: emission feature maximum channel=',
i7)
ier = 1
return
}
if (ifeattype == 1 & minch == 0) {
write (ttyout,304) minch
304 format ('ERROR: absorption feature minimum channel=',
i7)
ier = 1
return
}
#
# DEBUG:
#
# write (ttyout,*) 'DEBUG: bdmset: cl1=',cl1,' cl2=',cl2,
# ' cr1=',cr1,' cr2=',cr2
# do i = 1, cr2 {
# write (ttyout,*) 'DEBUG: bdmset: wav(',i,')=,wav(i),
# ' rflib(',i,')=,rflib(i),
# ' rflibc(',i,')=,rflibc(i)
# }
# write (ttyout,*) 'DEBUG: bdmset: minch=',minch,' ier=',ier
# write (ttyout,*) 'DEBUG: bdmset: maxch=',maxch
# END.DEBUG
#
# done!
#
return
end
U.S. Geological Survey,
a bureau of the U.S. Department of the Interior
This page URL= http://speclab.cr.usgs.gov/PAPERS/tricorder.1995/bdmset.r.html
This page is maintained by: Dr. Roger N. Clark rclark@speclab.cr.usgs.gov
Last modified November 18, 1998.