subroutine bandmp (wav,rflibc,rfobs,cl1,cl2,cr1,cr2,ictrol,minch,
maxch,ifeattype,iflag, ixpxl,iypxl,
rfobsc,xk,bd,rfit,slope,yintcp,rftemp,conref,
avlc,avrc)
implicit integer*4 (i-n)
#ccc name: bandmp
#ccc version date: January 29, 1990
#ccc author(s): Roger N. Clark
#ccc language: Ratfor
#ccc
#ccc short description:
#ccc This program does a least squares fit of a library
#ccc reference spectrum to another spectrum over a given wavelength
#ccc range.
#ccc
#ccc WARNING: Input setup must be correct (see input
#ccc parameters below). This is to maximize speed
#ccc in imaging spectroscopy applications.
#ccc
#ccc WARNING: subroutine bdmset MUST be called first!
#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 system requirements:
#ccc subroutines called:
#ccc
#ccc
#ccc parameter description:
#ccc INPUT:
#ccc wav = wavelengths R*4 (cr2 elements)
#ccc rflibc = reference library spectrum, continuum
#ccc removed R*4 (cr2 elements)
#ccc rfobs = observed 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 (not checked)
#ccc cr1 = continuum point begin on right side of band I*4
#ccc note: cr1 > cl2 + 1 (not checked)
#ccc cr2 = continuum point end on right side of band I*4
#ccc note: cr2 >= cr1 (not checked)
#ccc note: cr2 also determines the max array sizes
#ccc ictrol = error message control flag:
#ccc = 0 don't print error messages, just set
#ccc output to deleted values.
#ccc = 1 print error messages and set output to
#ccc deleted values.
#ccc minch = minimum channel in the library spectrum.
#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 iflag = flag to indicate if processing image-type data:
#ccc = 0 not image data (single spectrum analysis)
#ccc = 1 imaging spectrometer data.
#ccc ixpxl = x-pixel coordinate of image data (iflag = 1)
#ccc iypxl = y-pixel coordinate of image data (iflag = 1)
#ccc
#ccc OUTPUT:
#ccc rfobsc = continuum removed observed spectrum R*4 (cr2 elements)
#ccc NOTE: continuum is removed ONLY between cl2 and cr1
#ccc To get complete continuum, you must remove
#ccc it yourself using the slope and yintcp
#ccc variables below.
#ccc xk = k factor needed to make reference match library R*4
#ccc bd = band depth R*4
#ccc rfit = fit parameter normalized to 1 channel R*4
#ccc slope = slope to continuum of observed spectrum
#ccc yintcp = intercept to continuum of observed spectrum
#ccc rftemp = temporary working array R*4 (cr2 elements)
#ccc at a normal return, rftemp is the fitted,
#ccc continuum removed spectrum.
#ccc conref = continuum reflectance value of observed spectrum
#ccc avlc = left continuum reflectance value of observed spectrum
#ccc avrc = right continuum reflectance value of observed spectrum
#ccc
#ccc common description: none
#ccc message files referenced: none
#ccc internal variables:
#ccc suml = sum (rflibc) Real*8
#ccc sumll = sum (rflibc * rflibc) Real*8
#ccc sumol = sum (rfobsc * rflibc) Real*8
#ccc sumo = sum (rfobsc) Real*8
#ccc sumoo = sum (rfobsc * rfobsc) Real*8
#ccc file description: none
#ccc user command lines: none
#ccc update information:
#ccc NOTES:
#ccc
integer*4 cl1,cl2,cr1,cr2,ictrol,minch,maxch,ifeattype
real*4 wav(cr2),rflibc(cr2),rfobs(cr2)
real*4 rfobsc(cr2),xk,bd,rfit, rfit2
real*4 rftemp(cr2), slope, yintcp, xn
real*4 top, bottom, botm2, bprime, drfit2, drfit
real*4 avlc,avrc
real*8 suml,sumll,sumol,sumo,sumoo,dxk,dbb,dxn
real*8 dro, drl
integer*4 ttyout
ttyout = 6
#
# deleted point value:
#
delpt = -1.23e34
#
# first we must remove the continuum to the observed spectrum:
# compute rfobsc
#
# avlc = 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 || rflibc(i) == delpt ||
rfobs(i) == delpt) next
avlc = avlc + rfobs(i)
avwlc = avwlc + wav(i)
n = n +1
}
if (n < 1) {
if (ictrol == 1 & iflag == 0) write (ttyout,120)
120 format (' ERROR: all points in observed spectrum,',
' left continuum deleted.')
if (ictrol == 1 & iflag == 1) write (ttyout,121) ixpxl, iypxl
121 format (' ERROR: all points in observed spectrum,',
' left continuum deleted: pixel',i4,',',i4)
xk = delpt
bd = delpt
rfit = delpt
return
}
xn = float(n)
avlc = avlc / xn
avwlc = avwlc / xn
if (avlc < 0.1e-20) { # if the continuum id less than zero, delete output
if (ictrol == 1 & iflag == 0) write (ttyout,124)
124 format (' ERROR: left continuum <0.0.')
if (ictrol == 1 & iflag == 1) write (ttyout,125) ixpxl, iypxl
125 format (' ERROR: left continuum <0.0r: pixel',i4,',',i4)
xk = delpt
bd = delpt
rfit = delpt
return
}
# now right continuum average:
avrc = 0.0
avwrc = 0.0
n = 0
do i = cr1, cr2 {
if (wav(i) == delpt || rflibc(i) == delpt ||
rfobs(i) == delpt) next
avrc = avrc + rfobs(i)
avwrc = avwrc + wav(i)
n = n +1
}
if (n < 1) {
if (ictrol == 1 & iflag == 0) write (ttyout,122)
122 format (' ERROR: all points in observed spectrum,',
' right continuum deleted.')
if (ictrol == 1 & iflag == 1) write (ttyout,123) ixpxl, iypxl
123 format (' ERROR: all points in observed spectrum,',
' right continuum deleted: pixel',i4,',',i4)
xk = delpt
bd = delpt
rfit = delpt
return
}
xn = float(n)
avrc = avrc / xn
avwrc = avwrc / xn
if (avrc < 0.1e-20) { # if the continuum id less than zero, delete output
if (ictrol == 1 & iflag == 0) write (ttyout,126)
126 format (' ERROR: right continuum <0.0.')
if (ictrol == 1 & iflag == 1) write (ttyout,127) ixpxl, iypxl
127 format (' ERROR: right continuum <0.0r: pixel',i4,',',i4)
xk = delpt
bd = delpt
rfit = delpt
return
}
# now we have two pairs os x,y points to derive the continuum line
# thus: wav = a * rfobs + b
bottom = avwrc - avwlc
if (abs(bottom) < 0.1e-20) {
if (ictrol == 1) write (ttyout, 130)
130 format (' ERROR: wavelength range of continuum is',
' too small!')
xk = delpt
bd = delpt
rfit = delpt
return
}
a = (avrc - avlc)/bottom
b = avrc - a * avwrc
slope = a
yintcp = b
#
# now we can continuum correct the observed spectrum.
# compute rfobsc
#
# don't do because speed is reduced
#############
# if (cl1 > 1) {
# do i = 1, cl1-1 {
# rfobsc(i) = delpt
# }
# }
#############
# compute continuum only between cl2 and cr1 (but not including them)
# in imaging mode, but in single spectrum mode, include all
# of continuum.
if (iflag == 1) { # imaging mode, do minimum
il = cl2 + 1
ir = cr1 - 1
} else { # single spectrum mode, do more.
il = cl1
ir = cr2
}
do i = il, ir {
if (wav(i) == delpt || rflibc(i) == delpt ||
rfobs(i) == delpt) {
rfobsc(i) = delpt
} else {
contin = a * wav(i) + b
if (abs(contin) > 0.1e-20) {
rfobsc(i) = rfobs(i) / contin
} else {
rfobsc(i) = delpt
}
}
}
if (ifeattype == 1) {
if (wav(minch) == delpt) {
conref = delpt
} else {
conref = a * wav(minch) + b
}
} else {
if (wav(maxch) == delpt) {
conref = delpt
} else {
conref = a * wav(maxch) + b
}
}
#
# now compute sums between cl2 and cr1 (but not including them)
#
n = 0
suml = 0.0
sumll = 0.0
sumol = 0.0
sumo = 0.0
sumoo = 0.0
do i = il, ir {
if (rfobsc(i) == delpt) next
n = n +1
drl = dble(rflibc(i))
dro = dble(rfobsc(i))
suml = suml + drl
sumll = sumll + drl*drl
sumol = sumol + dro*drl
sumo = sumo + dro
# the following is needed only for correlation coefficient
sumoo = sumoo + dro*dro
}
if (n == 0) { # no sums, so delete output
if (ictrol == 1) write (ttyout, 140)
140 format (' ERROR: routine bandmp: sums in least ',
'squares have no data summed! Output deleted.')
xk = delpt
bd = delpt
rfit = delpt
return
}
#
# now compute a and b in the least squares solution
#
xn = float(n)
dxn = dble(xn)
top = sngl(sumol - (sumo * suml / dxn))
bottom = sngl(sumll - (suml * suml / dxn))
if (abs(bottom) < 0.1e-30) {
b = 0.0
} else {
b = top / bottom
}
#
# save the b value
#
bb = b
#
# now compute the xk value
#
if (abs(bb) < 0.1e-30) { # xk too large (b too small)
xk = delpt
bd = delpt
rfit = delpt
return
} else {
xk = (1.0 - b) / b
}
#
# now compute the reference spectrum that best matches the observed spectrum
#
xk1 = xk + 1.0
if (abs(xk1) < 0.1e-30) {
if (ictrol == 1 & iflag == 0) write (ttyout,222)
222 format (' ERROR: bandmp xk1 value too small,',
' output fit, depth deleted.')
if (ictrol == 1 & iflag == 1) write (ttyout,223) ixpxl, iypxl
223 format (' ERROR: bandmp xk1 value too small,',
' output fit, depth deleted: pixel',i4,',',i4)
bd = delpt
rfit = delpt
return
}
do i = il, ir {
if (rflibc(i) == delpt) {
rftemp(i) = delpt
} else {
rftemp(i) = (rflibc(i) + xk) / xk1
}
}
#
# compute band depth
#
if (ifeattype == 1) {
bd = 1.0 - rftemp(minch)
} else {
bd = 1.0 - rftemp(maxch)
}
#
# compute goodness of fit
#
# now compute aprime and bprime in the least squares solution
#
# top is still the same
botm2 = sngl(sumoo - (sumo * sumo / dxn))
if (abs(botm2) < 0.1e-30) {
bprime = 0.0
} else {
bprime = top / botm2
}
rfit2 = abs(bb * bprime)
rfit = sqrt(rfit2) # goodness of fit
#
# DEBUG:
#
# need to add: ictrol,ixpxl,iypxl,
#
# write (ttyout,*) 'DEBUG: bandmp: cl1=',cl1,' cl2=',cl2,
# ' cr1=',cr1,' cr2=',cr2
# do i = cl1, cr2 {
# write (ttyout,*) 'DEBUG: bandmp: wav(',i,')=',wav(i),
# ' rflibc(',i,')=',rflibc(i),
# ' rfobs(',i,')=',rfobs(i)
# }
# write (ttyout,*) 'DEBUG: bandmap: xn=',xn
# write (ttyout,*) 'DEBUG: bandmap: suml=',suml,
# ' sumll=',sumll
# write (ttyout,*) 'DEBUG: bandmap: sumol=',sumol,
# ' sumo=',sumo,
# ' sumoo=',sumoo
# write (ttyout,*) 'DEBUG: bandmap: xk=',xk,
# ' bd=',bd,' rfit2=',rfit2,' rfit=',rfit
# write (ttyout,*) 'DEBUG: bandmap: slope=',slope,
# ' yintcp=',yintcp
# write (ttyout,*) 'DEBUG: bandmap: top=',top,
# ' bottom=',bottom,
# ' botm2=',botm2
# write (ttyout,*) 'DEBUG: bandmap: bb=',bb,
# ' bprime=',bprime
# do i = cl1, cr2 {
# write (ttyout,*) 'DEBUG: bandmap: rfobsc(',i,')=',
# rfobsc(i), ' rftemp(',i,')=',rftemp(i)
# }
# write (ttyout,*) 'DEBUG: bandmp: minch=',minch,' iflag=',iflag
# END.DEBUG
#
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/bandmp.r.html
This page is maintained by: Dr. Roger N. Clark rclark@speclab.cr.usgs.gov
Last modified November 18, 1998.