#Original Function Location: /themis/lib/dav_lib/public_library/latest/spectral_tools.dvrc
/themis/lib/dav_lib/library/spectral_tools.dvrc
Source Code for Function: "make_band()" in "spectral_tools.dvrc" (Public)
spectral_tools_version=1.24
#1.01 version incremented for emcal rewrites
#1.02 added sma, sort_conc, sum_group_conc, summary_sma, and summary_sma_group
#1.03 minor edits to sma, sort_conc, sum_group_conc, summary_sma and summary_sma_group
#1.04 changes to sma (accepts a standard structure + resample option)
# changes to summary_sma, summary_sma_groups, now accept ranges (e.g. xlow=1, xhigh=10)
# dehyd2 changed to handle $DV_SCRIPT_FILES
#1.05 changes to sma to handle other data types than float
#1.06 added several generic spectral analysis tools
#1.07 added reform_speclib
#1.08 added i2i_tutorial and made make_band option to return a 73-pt spectrum
#1.09 dehyd modified to run on individual spectra in an XxY array
# sma - now returns elements that are all the same size as the endlib struct
# sma - modified exclude to take non-monotonically increasing values
# sma - fixed notchco2 to work on XxY arrays
#1.10 sma - added new elements to .grouped in sma output
#1.11 summary_sma(), summary_sma_group() - total rewrite, similar usage as before,
# summary_sma_group() now calls summary_sma() with a group=1 option
# modifictionas to sma() which allow for integration with new plot_sma and a bug fix
# to sum_group_conc() which now takes the precomputed covariance matrix from sma (fixes error values)
#1.12 sum_group_conc() - fixed goofy error where the number of elements was not being calculated correctly
#1.13 added tdb2struct - prc 11/2010
#1.14 sma - added nn option, which uses NNLS algorithm to fit spectra (faster and more accurate than original algorithm)
#1.15 fit_bb - added function to fit blackboides using nnls
#1.17 new versions of rad2tb2, tb2rad2, make_band - now take inst="instrument" options
# also included emissivity, cal_rad, irf_ri
#1.18 includes emcal2 and some major fixes to emssivity, cal_rad, and irf_ri
# also included r2t_swri
#1.19 fixed emcal2 to allow a list of environment temperatures, also data are now stored in x not y
#1.20 added sublib to allow the useer to subset all elements of a spectral library at once.
#1.21 fixed a bug in sma() that would not allow nn=1 to work with the surface=1 option.
#1.22 added structures to i2i() and it now works with MASTER as well.make_band() now has MASTER option as well
# modified emcal2() to have max_emiss option. Also fixed ctok() issues that were changed in emissivity
#1.23 added diffbb modeling option to fit_bb()
#1.24 bug fix to sma(); upgraded emcal2 to take new LabView spectrometer output
# norm_spec - taken from D. Rogers
# reform_speclib - created by c.edwards 8-5-10
# cal_rad - taken from /u/phil/scripts/dvinci/prc_beta.dv
# irf_ri - taken from /u/phil/scripts/dvinci/prc_beta.dv
# emissivity - taken from /u/phil/scripts/dvinci/prc_beta.dv
# tb2rad2 - taken from /u/phil/scripts/davinci/ir_analysis.dv
# rad2tb2 - taken from /u/phil/scripts/davinci/ir_analysis.dv
# make_band - taken from /u/phil/scripts/davinci/ir_analysis.dv
# make_temp_rad - taken from /u/phil/scripts/davinci/ir_analysis.dv
# i2i - taken from /u/phil/scripts/davinci/ir_analysis.dv
# i2i_tutorial - taken from /u/phil/scripts/davinci/ir_analysis.dv
# reverse_spectrum - taken from /u/phil/scripts/davinci/ir_analysis.dv
# cm2micron - taken from /u/phil/scripts/davinci/ir_analysis.dv
# read_vm - taken from /u/bandfield/josh.dvrc
# r2t_hi - taken from /u/baldridge/amb.dvrc
# r2t_lo - taken from /u/baldridge/amb.dvrc
# r2t_swri - taken from rt2_swri.dvrc - M.Osterloo
# emcal - taken from /u/baldridge/amb.dvrc
# emcal2 - created by c.edwards 3-17-11
# dehyd2 - taken from /u/cedwards/christopher.dvrc
# sublib - added by C. Edwards 4-8-2011
# sma - added by D. Rogers
# sort_conc - added by D. Rogers
# sum_group_conc - added by D. Rogers
# summary_sma - added by D. Rogers
# summary_sma_group - added by D. Rogers
# tdb2struct - added by P. Christensen
# fit_bb - added by c.edwards 1-24-2011
# sa - taken from /u/phil/scripts/atmos_rad.dv
define norm_spec(xaxis,w1,w2) {
if ($ARGC == 0) {
printf (" \n")
printf (" $1 = array of spectra to be normalized \n")
printf (" (xaxis = x-axis) default TES 73-point spectrum \n")
printf (" (w1 = wave1) default 200 cm-1 \n")
printf (" (w2 = wave2) default 1305 cm-1 \n")
printf (" \n")
return(0)
}
arry=$1
tesx=translate(xt2(),y,z)
axeis=cat(tesx[,,9:35],tesx[,,65:110],z)
wayve1=200
wayve2=1305
if(HasValue(xaxis)) {
axeis=xaxis
}
if(HasValue(w1)){
wayve1=w1
}
if(HasValue(w2)){
wayve2=w2
}
avearry=avg(arry,axis=xy)
dimarry=dim(arry)
nspecs=dimarry[1]
fspecs=create(dimarry[1],dimarry[2],dimarry[3],format=float,start=0,step=0)
for(i=1;i<=nspecs;i+=1){
strc=struct()
strc.data=arry[i]
strc.label="junk"
strc.xaxis=axeis
strc.shortlabel="shortjunk"
strc.id=0
fin1=sma(avearry,strc,wave1=wayve1,wave2=wayve2)
fspecs[i]=fin1.modeled
}
rtnstruc=struct()
rtnstruc.norm=fspecs
rtnstruc.ave=avg(fspecs,axis=xy)
rtnstruc.stddev=stddev(fspecs,axis=xy)
return(rtnstruc)
}
define reform_speclib(tes73,from,to) {
if($ARGC==0) {
printf("\n")
printf("Prepare a library from speclib.asu.edu for use with sma()\n")
printf("$1 = hdf downloaded from the spectral library tool\n\n")
printf("Optional:\n")
printf("tes = 1 crop to TES surface sensing channels only (73-point spectrum) (Default is 0)\n\n")
printf("from = direction to translate from (Default is \"y\")\n")
printf("to = direction to translate to (Default is \"x\")\n")
printf("c.edwards 8/5/2010\n\n")
return(null)
}
speclib=$1
if(dim(speclib.data)[1]==1) {
if(HasValue(from)==0) from="y"
if(HasValue(to)==0) to="x"
} else {
if(HasValue(from)==0) from="x"
if(HasValue(to)==0) to="x"
}
for(i=1;i<=length(speclib);i+=1) {
form=format(speclib[i])
if(form!="STRUCT" && form!="STRING" && form!="TEXT") {
speclib[i]=translate(speclib[i],from,to)
}
}
speclib.group=speclib.category
speclib.label=speclib.sample_name+" "+speclib.sample_id
speclib.shortlabel=speclib.sample_name
if(HasValue(tes73)==1) {
#crop to TES surface sensing channels only
speclib.xaxis=cat(speclib.xaxis[,,9:35],speclib.xaxis[,,65:110],axis=z)
speclib.data=cat(speclib.data[,,9:35],speclib.data[,,65:110],axis=z)
}
return(speclib)
}
define irf_ri(bb1, bb2, t1, t2, inst, xaxis, irf, noise_free, alice, riw1, riw2, instT, instdT) {
if(HasValue(bb1) ==0) {
printf(" Computes instrument response function (irf) and instrument radiance (ri) using\n")
printf(" two input blackbody spectra and temperatures. (Optionally uses 1 bb and an input irf)\n")
printf(" Usage: irf = irf_ri(bb1=spec1, t1=temp1, bb2=spec2, t2=temp2, inst='inst' [irf][xaxis=xaxis][noise_free=1][riw1][riw2][instT][instdT])\n")
printf(" t1 and t2 are in K\n")
printf(" Instruments supported: lab1, lab2, tes, tes5, tes73, mtes, themis, aster, tims \n")
printf(" [irf = spec] use with 1 blackbody (bb1 and t1) and a pre-existing irf\n")
printf(" [xaxis=xaxis] use xaxis instead of 'inst'\n")
printf(" Optional (for noise_free only)\n")
printf(" [noise_free = int] 0=no noise_free, 1=fit_bb noise_free, 2=fit diff bb noise_free (Default=0)\n")
printf(" [riw1 = f] wave1 for fit_bb range (default = 300)\n")
printf(" [riw2 = f] wave2 for fit_bb range (default = 300)\n")
printf(" [instT = f] set middle instrument temperature in fit_bb\n")
printf(" [instdT = f] set deviation around instrument temperature in fit_bb\n")
printf(" \n")
printf(" Ex: irf = irf_ri(bb1 = a.warmbb, bb2 = a.hotbb, t1 = ctok(a.warmbbt), t2 = ctok(a.hotbbt), inst='lab2')\n")
printf(" irf = irf_ri(bb1 = a.warmbb, bb2 = a.hotbb, t1 = ctok(a.warmbbt), t2 = ctok(a.hotbbt), inst='lab2', noise_free=1)\n")
printf(" irf = irf_ri(bb1 = a.warmbb, t1 = ctok(a.warmbbt), irf=irf_ri.irf, inst='lab2')\n")
printf("Returns .irf and .ri regardless of method used (e.g normal or noise_free)\n")
printf("p. christensen 1/20/2011\n")
return(0)
}
struc = struct()
if(HasValue(noise_free)==0) noise_free=0
# do initial tests for valid instrument or xaxis
if(HasValue(inst)) {
# user did specify an 'inst'
# test for valid values of inst
if(inst == "lab1" || inst == "lab2" || inst == "lab3" || inst == "tes" || inst == "tes5" || inst == "tes73" || inst == "mtes" || inst == "themis" || inst == "aster" || inst =="tims") {
# make the xaxis
xaxis = make_band(inst= inst)
struc.xaxis = xaxis
nbands = dim(xaxis)[3]
} else {
printf("not a valid instrument\n")
return(0)
}
} else {
# user did not specify an 'inst'
if(HasValue(xaxis)==0) {
printf("No instrument or xaxis specified\n")
return(0)
} else {
# no inst sent in, but do have an xaxis - assume it is a spectrometer
inst="spectrometer"
struc.xaxis = xaxis
}
}
struc.inst = inst
# do blackbody and temperature tests
if(HasValue(bb1)==0) {
printf("must send in at least 1 blackbody cal spectrum (and it should be in 'bb1'\n")
return(0)
} else {
lenbb1 = dim(bb1)[3]
if(HasValue(t1) == 0) {
printf("must send in at least 1 blackbody temperature (and it should be in t1)\n")
return(0)
} else {
# warn user to make sure temp is in kelvin
if(t1<200) {
printf("T1 seems low - make sure you input T in Kelvin\n")
#t1=ctok(t1)
}
struc.t1 = t1
# compute radiance of blackbody at temp t1
r1= bbr(xaxis, t1)
}
}
if(HasValue(bb2)==0) {
if(HasValue(irf) == 0) {
printf("if only send in 1 blackbody spectrum, must send in irf\n")
return(0)
}
} else {
lenbb2 = dim(bb2)[3]
if(lenbb2 != lenbb1) {
printf("length of bb1 spectrum not the same as bb2 spectrum\n")
return(0)
}
if(HasValue(t2)==0) {
printf("must send in second blackbody cal temperature\n")
return(0)
} else {
# warn user to make sure temp is in kelvin
if(t2<200) {
printf("T2 seems low - make sure you input T in Kelvin\n")
#t2 = ctok(t2)
}
struc.t2 = t2
# compute radiance of blackbody at temp t2
r2= bbr(xaxis, t2)
}
}
# test for correct sizes
lenxaxis = dim(xaxis)[3]
if(lenxaxis != lenbb1) {
printf("length of xaxis for this instrument and length of bb1 don't match\n")
return(0)
}
# compute instrument response function (irf)
if(HasValue(irf) == 0) {
irf = (bb2-bb1) / (r2 - r1)
}
# compute instrument radiance (ri)
ri = r1 - (bb1/irf)
if(noise_free == 0) {
# don't do the noise-free method
noise_free=0
struc.irf = irf
struc.ri = ri
# determine instrument temperture in default mode, even in not doing noise-free ri method
struc.noise_free_method = "none"
struc.inst_temp_method = "none"
} else if(noise_free==1 || noise_free==2){
struc.orig={}
struc.orig.ri=ri
struc.orig.irf=irf
# Compute noise free ri by fitting Planck function(s) to ri to get a noise-free ri
# and then compute a new irf using the noise-free ri with the average of the irf's
# computed from the two blackbody spectra
# the 'fit a library of bb's to ri' method
# this is the default (normal/appropriate) noise-free method
# set wavenumber range to fit over (300-2500 is default range)
if(HasValue(riw1)==0) {
riw1=300
}
if(HasValue(riw2)==0) {
riw2=2500
}
# fit a library of Planck functions to the radiance
if(noise_free==1) {
if(HasValue(instT)!=0 && HasValue(instdT)!=0) {
# use instT and instdT modes (see fit_bb documentation)
ri_struc = fit_bb(ri, xaxis, wave1=riw1, wave2=riw2, instT=instT, instdT=instdT)
struc.instT = instT
struc.instdT = instdT
} else {
# normal mode
ri_struc = fit_bb(ri, xaxis, type='inst', wave1=riw1, wave2=riw2)
}
} else if(noise_free==2) {
if(HasValue(instT)!=0 && HasValue(instdT)!=0) {
# use instT and instdT modes (see fit_bb documentation)
ri_struc = fit_bb(ri, xaxis, wave1=riw1, wave2=riw2, instT=instT, instdT=instdT, diffbb=1)
struc.instT = instT
struc.instdT = instdT
} else {
# normal mode
ri_struc = fit_bb(ri, xaxis, type='inst', wave1=riw1, wave2=riw2, diffbb=1)
}
}
nfri = ri_struc.modeled
# now go back and compute a noise-free irf
# use both blackbodies to determine irf,and then take the average
# this reduces the noise in the irf by sqrt(2)
irf1 = bb1 / (r1 - nfri)
irf2 = bb2 / (r2 - nfri)
nfirf=(irf1+irf2)/2.
# load structure
struc.irf = nfirf
struc.ri = nfri
struc.irf1 = irf1
struc.irf2 = irf2
if(noise_free==1) {
struc.noise_free_inst_temp = ri_struc.stats.avgtemp
struc.noise_free_method = "fit_bb"
} else if (noise_free==2) {
struc.noise_free_method = "fit_diffbb"
}
struc.riw1=riw1
struc.riw2=riw2
} else if (noise_free==3) {
# This version is for testing only to compare with Alice's version used in emcal
# Do not use this option !!!!
# Mimic Alice Baldridge's approach used in emcal
# do a real simple version - find ave_t from middle 1/3 of spectrum and use that
# to compute the radiance from the Planck function
bt = rad2tb2(ri, xaxis=xaxis)
nz=dim(bt)[3]
z1=int(nz*.1)
z2=int(nz*.26)
avet=avg(bt[,,z1:z2])
nfri = tb2rad2(avet, xaxis=xaxis)
struc.ri = nfri
struc.irf = irf
struc.ave_t_w1 = xaxis[,,z1]
struc.ave_t_w2 = xaxis[,,z2]
struc.noise_free_inst_temp = avet
struc.noise_free_method = "alice"
}
return(struc)
}
define cal_rad(spec, irf, ri, inst, xaxis) {
if(HasValue(spec) == 0) {
printf(" Compute calibrated radiance\n")
printf(" Usage: rad = cal_rad(spec=radiance_spectrum, irf=irf, ri=ri, inst='inst' [or xaxis=xaxis])\n")
printf(" Use inst = 'inst' mode if standard instrument, or send in xaxis\n")
printf(" Instruments supported: lab1, lab2, tes, tes5, tes73, mtes, themis, aster, tims\n")
printf(" Ex: rad = cal_rad(spec = a.data, irf = irf_ri.irf, ri = irf_ri.ri, inst='lab2')\n")
printf(" rad = cal_rad(spec = a.data, irf = irf_ri.irf, ri = irf_ri.ri, xaxis=a.xaxis)\n")
printf(" p. christensen 1/20/2011\n")
return(0)
}
struc = struct()
if(HasValue(spec) == 0) {
printf("no spectra sent in\n")
return(0)
}
if(HasValue(irf) == 0) {
printf("no instrument response function (irf) sent in\n")
return(0)
}
if(HasValue(ri) == 0) {
printf("no instrument radiance (ri) sent in\n")
return(0)
}
if(HasValue(inst) == 0) {
# user did not specify an 'inst' - must then provide an xaxis
if(HasValue(xaxis)==0) {
printf("No instrument or xaxis specified\n")
return(0)
} else {
# no inst sent in, but do have an xaxis - assume it is a spectrometer
inst = "spectrometer"
struc.xaxis = xaxis
}
} else {
# user did specify an 'inst' - test for valid values of inst
if(inst == "lab1" || inst == "lab2" || inst == "lab3" || inst == "tes" || inst == "tes5" || inst == "tes73" || inst == "mtes" || inst == "themis" || inst == "aster" || inst =="tims") {
# make the xaxis
xaxis = make_band(inst= inst)
struc.xaxis = xaxis
nbands = dim(xaxis)[3]
} else {
printf("not a valid instrument\n")
return(0)
}
}
struc.inst = inst
# test for correct sizes
lenirf = dim(irf)[3]
lenri = dim(ri)[3]
lenxaxis = dim(xaxis)[3]
if(lenxaxis != lenirf || lenxaxis != lenri || lenirf != lenri) {
printf("length of xaxis for this instrument, length of irf, and/or length of ri don't match\n")
return(0)
}
# compute calibrated radiance
radiance = (spec/irf) + ri
struc.radiance = radiance
return(struc)
}
define emissivity(inst, xaxis, wave1, wave2, filter, b1, b2, max_emiss, max_emiss_low, downwelling_e, downwelling_t, downwelling_rad, dehyd, toffset) {
# compute emissivity from calibrated spectral radiance
# resurrected - prc 12/15/2010 during Fall AGU. Name changed to 'emissivity3' (Changed back to emissivity2 1/30/2011)'
# modified from themis_emissivity prc 12/5/08
# p. christensen 2/18/02
# returns input structure with emissivity and maximum brightness temperature
# normally run this function after geometric registration has been done to align bands
if ($ARGC == 0) {
printf (" \n")
printf (" Generic script to compute emissivity from standard calibrated IR spectral radiance datasets\n")
printf (" \n")
printf (" usage: e = emissivity(radiance, inst=[lab1, lab2, tes, tes5, mtes, themis, aster, tims], [band_list]\n")
printf (" [wave1][wave2][b1][b2][filter]\n")
printf (" [max_emiss][max_emiss_low][downwelling_e][downwelling_t][downwelling_rad])\n")
printf (" returns structure with emissivity (emissivity), target_temperature (maxbtemp), and parameters used\n")
printf (" Instruments supported: lab1, lab2, tes, tes5, mtes, themis, aster, tims \n")
printf (" \n")
printf (" [downwelling_t=f] temperature of downwelling radiance (e.g. environmental chamber)\n")
printf (" [downwelling_e=f] emissivity of downwelling radiance (e.g. environmental chamber)(default=1.)\n")
printf (" [downwelling_rad=spec] spectrum of downwelling radiance\n")
printf (" [dehyd = 1] use dehy2 function\n")
printf (" [filter = n] filter size for determining max brightness temperature \n")
printf (" [max_emiss = f] optional maximum emissivity (e.g. 0.98) (default = 1.) \n")
printf (" [max_emiss_low = f] optional maximum emissivity used for low-T calc (for TES) (default = 1.) \n")
printf (" THEIMS: band_list is required (n,1,1) array of the n band numbers contained in the input radiance array)\n")
printf (" spectrometers: [wave1 = n] optional start wavenumber for brightness temperature calc (each has default) \n")
printf (" spectrometers: [wave2 = n] optional end wavenumber for brightness temperature (each has default)\n")
printf (" broadband: [b1 = n] optional start band for brightness temperature calculation\n")
printf (" broadband: [b2 = m] optional end band for brightness temperature calculation\n")
printf (" \n")
printf (" THEMIS: Default uses Band 3 through Band 9 to determine maximum brightness temp for 10-band image\n")
printf (" THEMIS: Default uses image planes 1-2 to determine maximum brightness temp for 3-band image (assumes Bands 4, 9, 10)\n")
printf (" \n")
printf (" \n")
printf (" examples: e = emissivity(rad, inst='mtes')\n")
printf (" e = emissivity(themis=1, a.spectral_qube.data, a.spectral_qube.band_bin.band_bin_band_number)\n")
printf (" (where a is RDR structure loaded using davinci 'load_pds' function)\n")
printf (" \n")
printf (" e = emissivity(themis=1, a.spectral_qube.data, a.spectral_qube.band_bin.band_bin_band_number)\n")
printf (" e = emissivity(themis=1, a.spectral_qube.data, a.spectral_qube.band_bin.band_bin_band_number, max_emiss = 0.98)\n")
printf (" e = emissivity(themis=1, rad, 1//2//3//4//5//6//7//8//9//10, b1 = 4, b2 = 8)\n")
printf (" e = emissivity(themis=1, rad, band_list)\n")
printf (" \n")
printf (" e = emissivity(a.radiance, inst='lab2', downwelling_t = ctok(a.envtemp))\n")
printf (" \n")
printf (" Notes: For THEMIS and other imagers should only be done after geometric registration. \n")
printf (" Function now sets radiance values <0. to 0. using: rad[where(rad<0.)] = 0. \n")
printf (" \n")
printf (" Generic version: p. christensen 12/08\n")
printf (" THEMIS version: p. christensen 4/02-4/03\n")
printf (" Mini-TES version: p. christensen 1/04-9/04\n")
return(0)
}
rad = $1
size=dim(rad)
struc = struct()
# check for xaxis, rather than inst
if(HasValue(inst) == 0) {
# no inst sent in
if(HasValue(xaxis)==0) {
printf("No inst or xaxis sent in - sorry\n")
return(0)
}
inst = 'spectrometer'
xaxis = xaxis
}
######## Step 1 - Step up parameters *************************
# check for spectrometer or broadband
if(inst=='lab1' || inst=='lab2' || inst=='lab3' || inst=='tes' || inst=='tes5' || inst=='mtes' || inst=='spectrometer') {
# spectrometers
if(inst=='lab1' || inst=='lab2' || inst=='lab3') {
# ASU lab spectrometers
# make xaxis
xaxis = make_band(inst=inst)
nbands = dim(xaxis)[3]
# set parameters for target temperature determination
# set default wavenumber range over which to find maximum brightness temperature
if(HasValue(wave1) == 0) {
wave1 = 500
}
if(HasValue(wave2) == 0) {
wave2 = 1700
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
max_emiss = 1.
}
if(HasValue(max_emiss_low) == 0) {
max_emiss_low = max_emiss
}
# threshold temperature and wavenumber range for
# low vs full wavenumber range to determine brightness temp
threshold_t_cold = 160.
threshold_t_warm = 170.
w1_cold = 450
w2_cold = 900
# set size of filter used to filter brightness temperature
if(HasValue(filter) == 0) {
filter = 3
}
# set up for CO2-ignore region
co2chan = 0
co2lowchan = 0
co2highchan = 0
# set up downwelling radiance
if(HasValue(downwelling_t) == 0) {
downwelling_t = 0.
}
if(HasValue(downwelling_e) == 0) {
downwelling_e = 1.
}
if(HasValue(downwelling_rad) == 0) {
downwelling_rad = 0.
}
} else if(inst=='tes') {
# from TES Data Processing User's Guide - 2006 version
# set parameters for target temperature determination
# set default wavenumber range over which to find maximum brightness temperature
xaxis = make_band(inst = "tes")
nbands = dim(xaxis)[3]
if(HasValue(wave1) == 0) {
wave1 = 300
}
if(HasValue(wave2) == 0) {
wave2 = 1350
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
max_emiss = 1.
}
if(HasValue(max_emiss_low) == 0) {
# TES processing sets max emissivity for low-T to 0,97
max_emiss_low = 0.97
}
# threshold temperature and wavenumber range for
# low vs full wavenumber range to determine brightness temp
threshold_t_cold = 215.
threshold_t_warm = 225.
w1_cold = 300
w2_cold = 500
# set size of filter used to filter brightness temperature
if(HasValue(filter) == 0) {
filter = 7
}
# set up for CO2-ignore region
# TES calibration excluded the region from '500 to 800 cm-1' (TES Data
# Processing User's Guide)
# modified from Mini-TES version to give correct index
# set up to be 500 to 800 cm-1
co2chan = maxchan(int(xaxis - 667) == 0)
co2lowchan = co2chan - 15
co2highchan = co2chan + 12
# set up downwelling radiance
downwelling_t = 0.
downwelling_e = 0.
downwelling_rad = 0.
} else if(inst=='tes5') {
# set parameters for target temperature determination
# set default wavenumber range over which to find maximum brightness temperature
xaxis = make_band(inst="tes5")
nbands = dim(xaxis)[3]
if(HasValue(wave1) == 0) {
wave1 = 300
}
if(HasValue(wave2) == 0) {
wave2 = 1350
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
max_emiss = 1.
}
if(HasValue(max_emiss_low) == 0) {
# TES processing sets max emissivity for low-T to 0,97
max_emiss_low = 0.97
}
# threshold temperature and wavenumber range for
# low vs full wavenumber range to determine brightness temp
threshold_t_cold = 215.
threshold_t_warm = 225.
w1_cold = 300
w2_cold = 500
# set size of filter used to filter brightness temperature
if(HasValue(filter) == 0) {
filter = 7
}
# set up for CO2-ignore region
# TES calibration excluded the region from '500 to 800 cm-1' (TES Data
# Processing User's Guide)
# noel's version
# modified from Mini-TES version to give correct index
# set up to be 500 to 800 cm-1
co2chan = maxchan(int(xaxis - 667) == 0)
co2lowchan = co2chan - 31
co2highchan = co2chan + 25
# set up downwelling radiance
downwelling_t = 0.
downwelling_e = 0.
downwelling_rad = 0.
} else if(inst=='mtes') {
# set parameters for target temperature determination
# set default wavenumber range over which to find maximum brightness temperature
xaxis = make_band(inst = "mtes")
nbands = dim(xaxis)[3]
if(HasValue(wave1) == 0) {
wave1 = 500
}
if(HasValue(wave2) == 0) {
wave2 = 1400
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
max_emiss = 1.
}
if(HasValue(max_emiss_low) == 0) {
max_emiss_low = max_emiss
}
# threshold temperature and wavenumber range for
# low vs full wavenumber range to determine brightness temp
threshold_t_cold = 220.
threshold_t_warm = 230.
w1_cold = 450
w2_cold = 900
# set size of filter used to filter brightness temperature
if(HasValue(filter) == 0) {
filter = 7
}
# set up for CO2-ignore region
# noel's version
# this is what is used in PRC's MER Mini-TES calibration
# modified from Mini-TES version to give correct index
# set up to be 500 to 800 cm-1
co2chan = maxchan(int(xaxis - 669) == 0)
co2lowchan = co2chan - 10
co2highchan = co2chan + 10
# set up downwelling radiance
downwelling_t = 0.
downwelling_e = 0.
downwelling_rad = 0.
} else if(inst=='spectrometer') {
# use lab2 default values
if(HasValue(wave1) == 0) {
wave1 = 500
}
if(HasValue(wave2) == 0) {
wave2 = 1700
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
max_emiss = 1.
}
if(HasValue(max_emiss_low) == 0) {
max_emiss_low = max_emiss
}
# threshold temperature and wavenumber range for
# low vs full wavenumber range to determine brightness temp
threshold_t_cold = 160.
threshold_t_warm = 170.
w1_cold = 450
w2_cold = 900
# set size of filter used to filter brightness temperature
if(HasValue(filter) == 0) {
filter = 3
}
# set up for CO2-ignore region
co2chan = 0
co2lowchan = 0
co2highchan = 0
# set up downwelling radiance
if(HasValue(downwelling_t) == 0) {
downwelling_t = 0.
}
if(HasValue(downwelling_e) == 0) {
downwelling_e = 1.
}
if(HasValue(downwelling_rad) == 0) {
downwelling_rad = 0.
}
}
######## Step 2 - Determine target temperature ##################
struc.xaxis=xaxis
# compute brightness temperature
# scale radiance to maximum emissivity (default = 1.)
tb = btemp(xaxis, rad*1./max_emiss)
# compute target (kinetic) temperature
# copy brightness temperture to work array
t1 = tb
# mask out CO2 region if desired
if(co2lowchan == 0 && co2highchan == 0) {
# do nothing if neither CO2 range set
} else {
# zero out CO2 region for finding TB
t1[,,co2lowchan:co2highchan] = 0.
}
# find maximum brightness temperatue
# first step - use full wavenumber range
# set values outside wavenumber range used to find maximum brightness temp to 0.
t1[where(xaxis < wave1)] = 0.
t1[where(xaxis > wave2)] = 0.
# filter brightness temperature to reduce noise
kernel = create(1,1,filter,start=1, step=0)
tfilt1 = convolve(t1, kernel, 1)
#tfilt1 = boxfilter(t1, filter)
# find maximum brightness temperature
max_t1 = max(tfilt1, axis=z)
nx = dim(tfilt1)[1]
ny = dim(tfilt1)[2]
pos_t1 = create(nx, ny, 1, format=int)
for(i=1; i<=nx; i+=1) {
for(j=1; j<=ny; j+=1) {
pos_t1[i,j] = maxpos(tfilt1[i,j])[3]
}
}
# second step - use cold-temperature wavenumber range
if(inst=='tes' || inst=='tes5') {
if(max_emiss == max_emiss_low) {
# don't do unnecessary work if emissivity values are the same
t2 = tb
} else {
t2 = btemp(xaxis, rad*1./max_emiss_low)
}
} else {
# this was the algorithm for mini-TES
t2 = tb
}
# mask out CO2 region if desired
if(co2lowchan == 0 && co2highchan == 0) {
# do nothing if neither CO2 range set
} else {
# zero out CO2 region for finding TB
#### this step was not done in original mtes script (emissivity()) prc 2/5/09
t2[,,co2lowchan:co2highchan] = 0.
}
# mask out wavenumber region exclude for cold case
t2[where(xaxis < w1_cold)] = 0.
t2[where(xaxis > w2_cold)] = 0.
# filter brightness temperature to reduce noise
tfilt2 = convolve(t2, kernel, 1)
#tfilt2 = boxfilter(t2, filter)
# find maximum brightness temperatue
max_t2 = max(tfilt2, axis=z)
pos_t2 = create(nx, ny, 1, format=int)
for(i=1; i<=nx; i+=1) {
for(j=1; j<=ny; j+=1) {
pos_t2[i,j] = maxpos(tfilt2[i,j])[3]
}
}
# loop through array - check each point separately
maxbtemp = max_t1 * 0.
nx = dim(maxbtemp)[1]
ny = dim(maxbtemp)[2]
for(i=1; i<=nx; i+=1) {
for(j=1; j<=ny; j+=1) {
if(max_t1[i,j] > threshold_t_warm) {
# this was a warm target
maxbtemp[i,j] = max_t1[i,j]
} else if (max_t2[i,j] < threshold_t_cold) {
# this was a cold target
maxbtemp[i,j] = max_t2[i,j]
} else {
# make smooth transition around threshold_t
maxbtemp[i,j] = (max_t1[i,j] + max_t2[i,j]) / 2.
}
}
}
# compute radiance of maximum brightness temperature
rad_bb = rad * 0.
for(j=1; j<= nx; j+=1) {
for(k=1; k<= ny; k+=1) {
if(HasValue(toffset)) {
printf("offseting target temperature by %.2f\n", toffset)
bt = maxbtemp[j,k,]+toffset
} else {
bt = maxbtemp[j,k,]
}
rad_bb[j,k,] = bbr(xaxis, bt)
#rad_bb[j,k,] = bbr(xaxis, maxbtemp[j,k,])
}
}
######## compute emissivity ##################
# include downwelling radiance if desired
if(downwelling_t == 0 && downwelling_rad == 0) {
# no downwelling radiance
printf("no downwelling radiance\n")
emiss = rad / rad_bb
} else {
# there is downwelling radiance
printf("doing downwelling radiance\n")
# compute downwelling radiance
if(downwelling_rad == 0.) {
downwelling_rad = downwelling_e * bbr(clone(xaxis,x=dim(downwelling_t)[1],y=dim(downwelling_t)[2]), downwelling_t)
}
emiss = (rad - downwelling_rad) / (rad_bb - downwelling_rad)
}
######## remove atmospheric water bands (mostly for lab measurements)
if(HasValue(dehyd)) {
emiss = dehyd2(emiss, xaxis)
struc.dehyd = dehyd
} else {
struc.dehyd = 0
}
# set bad values to 0.
emiss[where(rad<=0.)]=0.
# add items to structure that are unique to spectrometer
# keep mtes naming convention
struc.wave1 = wave1
struc.wave2 = wave2
struc.w1_cold = w1_cold
struc.w2_cold = w2_cold
struc.threshold_t_cold = threshold_t_cold
struc.threshold_t_warm = threshold_t_warm
struc.filter = filter
struc.co2chan = co2chan
struc.co2lowchan = co2lowchan
struc.co2highchan = co2highchan
struc.tb = tb
struc.tfilt1 = tfilt1
struc.tfilt2 = tfilt2
struc.max_t1 = max_t1
struc.max_t2 = max_t2
struc.pos_t1 = pos_t1
struc.pos_t2 = pos_t2
struc.rad_bb = rad_bb
struc.max_emiss_low = max_emiss_low
} else {
######## do broadband instruments ####################
# broadband instrument
# none of the broadband instruments compute downwelling radiance correction
downwelling_t = 0.
downwelling_e = 0.
downwelling_rad = 0.
if(inst=='aster') {
# assume all 5 bands are sent
if(size[3] != 5) {
printf("Not all 5 ASTER bands sent in - sorry\n")
return(0)
}
band_list = 1//2//3//4//5
if(HasValue(b1)==0) {
b1=1
}
if(HasValue(b2)==0) {
b2 = 5
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
# set to standard value used by ASTER team (prc - check 1/6/09)
max_emiss = 0.97
}
} else if(inst=='tims') {
# assume all 6 bands are sent
if(size[3] != 6) {
printf("Not all 6 TIMS bands sent in - sorry\n")
return(0)
}
band_list = 1//2//3//4//5//6
if(HasValue(b1)==0) {
b1=1
}
if(HasValue(b2)==0) {
b2 = 6
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
# set to standard value used by ASTER team (prc - check 1/6/09)
max_emiss = 0.97
}
} else if(inst=='modis') {
} else {
# default is themis
if($ARGC == 1) {
# no band list - check to make sure correct number of bands
if(size[3] != 10 && size[3] != 3) {
printf("Problem - incorrect number of bands (%d) without a band list\n", size[3])
return(0)
} else {
if(size[3] == 10) {
# sent in full image
band_list = 1//2//3//4//5//6//7//8//9//10
# set bands to use to determine max brightness temperature
if(HasValue(b1)==0) {
b1=3
}
if(HasValue(b2)==0) {
b2 = 9
}
} else {
# set up for normal 3-band nighttime image
band_list = 3//8//9
# set bands to use to determine max brightness temperature
if(HasValue(b1)==0) {
b1 = 1
}
if(HasValue(b2)==0) {
b2 = 2
}
}
}
} else if($ARGC == 2) {
# sent in a band list
band_list = $2
# set bands to use to determine max brightness temperature
if(HasValue(b1)==0) {
b1 = 1
}
if(HasValue(b2)==0) {
b2 = dim($1)[3]
}
} else {
# a problem
printf("Problem - incorrect number of arguments (%d)\n", $ARGC)
}
# set maximum emissivity values
if(HasValue(max_emiss) == 0) {
max_emiss = 1.
}
}
# set rad values < 0. to 0.
rad[where(rad<0.)] = 0.
# compute brightness temperature - adjust measured radiance assuming to corresponds to a
# maximum emissivity (default = 1.)
btemp = rad2tb2( (rad/max_emiss), band_list, inst=inst)
btemp[where(btemp<100.)] = 0.
# compute emissivity of scene
# find maximum brightness temperature - assume to be kinetic temperature
printf("\nfinding max brightness temperature using bands %d to %d, maximum emissivity = %.2f\n\n", b1, b2,max_emiss)
maxbtemp = max(btemp[,,b1:b2], axis = z)
emiss = btemp * 0.
rad_max_btemp = btemp * 0.
# compute radiance of the maximum brightness temp for all bands
rad_max_btemp = tb2rad2(maxbtemp, inst=inst)
# compute emissivity - divide by greybody
emiss = rad / rad_max_btemp
struc.band_list = band_list
struc.band1 = b1
struc.band2 = b2
struc.btemp = btemp
}
struc.downwelling_e = downwelling_e
struc.downwelling_t = downwelling_t
struc.downwelling_rad = downwelling_rad
struc.max_emiss = max_emiss
struc.emissivity = emiss
struc.target_temperature = maxbtemp
#struc.maxbtemp = maxbtemp
return(struc)
}
define rad2tb2(inst, max_emiss, lab1, lab2, lab3, tes, tes5, tes73, mtes, themis, aster, tims, modis, xaxis) {
if ($ARGC == 0) {
printf (" \n")
printf (" General routine to convert spectral radiance to temperature \n")
printf (" For spectrometers units the are W cm-2 sr-1 /cm-1\n")
printf (" For broadband radiometers the units are W cm-2 sr-1 micron-1) \n")
printf (" Temperature array returned is in same order as wavelength or wavenumber\n")
printf (" of instrument (e.g. ascending wavenumber or wavelength) \n")
printf (" Default instrument is THEMIS\n")
printf (" \n")
printf (" usage: temp = rad2tb2(radiance, inst=[lab1,lab2, lab3, tes, tes5, mtes, themis, aster, tims, modis] [max_emiss])\n")
printf (" Input 1 radiance and optional instrument (default is THEMIS)\n")
printf (" or: temp = rad2tb2(rad_array, inst='tes', 'themis', 'aster', etc)\n")
printf (" Input one radiance and instrument\n")
printf (" Returns temp[1,1,nbands]\n")
printf (" \n")
printf (" or: temp = rad2tb2(rad_array, inst=themis, aster, etc))\n")
printf (" Input array of radiances in nbands: rad_array[nx, ny, nbands] \n")
printf (" Dimension of nbands must match instrument \n")
printf (" Returns temp[nx, ny, nbands]\n")
printf (" \n")
printf (" or: temp = rad2tb2(rad_array, bandlist, inst='themis', 'aster', etc]))\n")
printf (" Input array of radiances in partial set of 1 or more bands: rad_array[nx, ny, ibands] \n")
printf (" Broadband instruments only\n")
printf (" Dimension of ibands must match # of bands in bandlist \n")
printf (" Example: temp = rad2tb2(image, cat(3,5,7,x), inst=themis) \n")
printf (" Example: temp = rad2tb2(cat(.00041, .00042, .00043,z), cat(4,5,6,x), inst='themis') \n")
printf (" Returns temp[nx, ny, ibands]\n")
printf (" \n")
printf (" [max_emiss = f] Set emissivity to value other than 1.0 (default)\n")
printf (" \n")
printf (" Note: Reads in appropriate radiance-to-temperature look-up table\n")
printf (" \n")
return(0)
}
# modified to include both 'tes=1' and 'inst="tes"' modes. prc `1/28/2011.
# changed from rad2tb2 to inst='lab1' mode, prc 12/16/2010
# For spectrometers this routine accepts:
# 1 radiance or an array of rads (z dim must equal # bands (e.g. 143 for TES))
# For broadband instrument this routine accepts:
# 1 radiance, an array of radiances, an array of radiances in all bands for the instrument,
# or an array of radiances in a subset of bands with a list of bands
# p. christenen 11/15/08 - derived from tb2rad2
# and tbtorad - p. christensen 10/2/01
if(HasValue(xaxis)) {
inst = 'spectrometer'
} else {
# check for 'old way' (tes=1) and convert to 'new way' (inst="tes")
args=cat("lab1", "lab2", "lab3", "tes", "tes5", "tes73", "mtes", "themis", "aster", "tims", "modis", axis=y)
for(i=1; i<=length(args); i+=1) {
if(HasValue(eval(args[,i]))==1) {
inst = args[,i]
}
}
}
# check for spectrometer or broadband
if(inst=='lab1' || inst=='lab2' || inst=='lab3' || inst=='tes' || inst=='tes5' || inst=='mtes' || inst=='spectrometer') {
# spectrometer
rad = $1
if(HasValue(xaxis)) {
x = xaxis
} else {
x = make_band(inst=inst)
}
nbands = dim(x)[3]
size=dim(rad)
temp = create(size[1], size[2], nbands, format=float)
# loop over x, y elements in input rad array
# function brtemp checks for temps < 0
for(i=1; i<= size[1]; i+=1) {
for(j=1; j<= size[2]; j+=1) {
if(size[3] == 1) {
# sent in 1 temp
temp[i,j] = btemp(x, rad[i,j,1])
} else {
# sent in array of rads - 1 for each band
# z dimension must match nbands
if(size[3] != nbands) {
printf("z dimension (%d) of input rad array does not match nbands (%d)\n", size[3], nbands)
return(0)
} else {
temp[i,j] = btemp(x, rad[i,j])
}
}
}
}
return(temp)
}
# read in temp_to_radiance (tem_rad)array made by 'make_temp_rad' for broadband instruments
# temp_rad arrays are dimensioned: [xaxis+nbands, ntemps, 1]
# they contain the wavelength(microns) and nbands of integrated radiance values for each temperature
# default instrument is themis
# default path set up by Chris Edwards for davinci use on multiple platforms
PATH = $DV_SCRIPT_FILES
PATHI = PATH + "/instrument_parameters"
if(inst=='aster') {
temp_rad =read(PATHI+"/temp_rad_aster.vic")
} else if(inst=='tims') {
temp_rad =read(PATHI+"/temp_rad_tims.vic")
} else if(inst=='modis') {
# not done yet
} else if(inst=='themis') {
temp_rad = read(PATHI+'/temp_rad_v4')
} else {
# default is themis
temp_rad = read(PATHI+'/temp_rad_v4')
}
# set up the broadband instruments
nbands = dim(temp_rad)[1] - 1
ntemps = dim(temp_rad)[2]
start_temp = temp_rad[1,1]
end_temp = temp_rad[1,ntemps]
if($ARGC == 1) {
# either 1 rad, an array of radiances, or
# an nx x ny array in nbands, where nbands must match instrument
# no band list
rad = $1
size = dim(rad)
if(size[3] == 1) {
# 1 rad or an array of rads
if(size[1]== 1 && size[2] == 1) {
# correct size for 1 radiance
# clone temperature to match number of bands in temp_rad array
rad_array = clone(rad, 1, 1, nbands)
} else if(size[1] != 1 || size[2] != 1) {
# an array of rads with only 1 band and no list specified
rad_array = clone(rad,1, 1, nbands)
}
} else {
# sent in an array of temps that must be dimensioned [nx, ny, nbands]
if(size[3] != nbands) {
# a problem
printf("incorrect z-size - must match number in instrument\n")
return(0)
} else {
#printf("ok - %d by %d array with correct number of bands %d\n", size[1], size[2], size[3])
rad_array = rad
}
}
# using all the bands in the instrument - set temp_rad_interp to temp_rad
temp_rad_interp = temp_rad
} else if ($ARGC == 2) {
# sent in an array of temperatures in 1 or more bands and specified which band(s)
rad_array = $1
size = dim(rad_array)
bandlist = $2
list_size = dim(bandlist)
# check format of bandlist
if(list_size[2] >1 || list_size[3] > 1) {
printf("bad format for bandlist - bands should be in x; ny and nz should be 1\n")
return(0)
}
nlist = list_size[1]
if(size[3] != nlist) {
# a problem
printf("incorrect z-size - if specifing 1 or more bands, nbands in image must be same as nbands in list\n")
return(0)
}
if(nlist > nbands) {
# a problem
printf("incorrect band choice - higher than number of bands in this instrument\n")
return(0)
}
# test to make sure no band in the list is higher than # of bands
max_list = max(bandlist[,1,1])
if(max_list > nbands) {
# a problem
printf("a band in the list (%d) is higher than the number of bands (%d)\n", max_list, nbands)
return(0)
}
# ok - passed the tests
if(nlist == 1) {
# only using 1 band: move the correct band from temp_rad array into temp_rad_interp
temp_rad_interp = cat(temp_rad[1], temp_rad[1+bandlist], x)
} else {
# sent in a list of bands - move the appropriate bands of temp_rad into temp_rad_interp
temp_rad_interp = create(nlist+1, ntemps, 1, format=float)
temp_rad_interp[1] = temp_rad[1]
for(i=1; i<=nlist; i+=1) {
temp_rad_interp[i+1] = temp_rad[bandlist[i] + 1]
}
}
}
# now do the easy part - interpolate in radiance vs temperature array
nbands_used = dim(temp_rad_interp)[1]-1
# debug mode
#printf("nbands=%d ntemps=%d size[1]=%d size[2]=%d size[3]=%d nbands_used=%d\n", nbands, ntemps, size[1], size[2], size[3], nbands_used)
temp = rad_array * 0.
# scale radiance to max emissivity
if(HasValue(max_emiss) == 0) {
max_emiss = 1.
}
#printf("max emiss = %.2f\n", max_emiss)
rad_array = rad_array / max_emiss
for(i=1; i <= nbands_used; i+=1) {
index = i+1
temp[,,i] = interp(temp_rad_interp[1], temp_rad_interp[index], rad_array[,,i])
}
# clean up to free memory
rad_array = 0.
# final test for temperatures outside range in temp_rad arrays
temp[where(temp < start_temp)] = 0.
temp[where(temp > end_temp)] = 0.
return(temp)
}
define tb2rad2(inst, lab1, lab2, lab3, tes, tes5, tes73, mtes, themis, aster, tims, modis, xaxis) {
if ($ARGC == 0) {
printf (" \n")
printf (" General routine to convert temperature to spectral radiance \n")
printf (" For spectrometers the units the are W cm-2 sr-1 /cm-1\n")
printf (" For broadband radiometers the units are W cm-2 sr-1 micron-1\n")
printf (" Radiance array returned is in same order as input wavelength or wavenumber\n")
printf (" of instrument (e.g. ascending wavenumber or wavelength) \n")
printf (" Default instrument is THEMIS\n")
printf (" \n")
printf (" usage: rad = tb2rad2(temp(K), inst=[lab1, lab2, tes, tes5, mtes, themis, aster, tims, modis])\n")
printf (" Input one temperature and instrument (default is THEMIS)\n")
printf (" inst='lab1' etc - specifies which instrument\n")
printf (" Returns rad[1,1,nbands]\n")
printf (" \n")
printf (" or: rad = tb2rad2(temp_array(K), inst='lab1',etc)\n")
printf (" Input array of temps and instrument: temp_array[nx, ny]\n")
printf (" Returns rad[nx,ny,nbands] \n")
printf (" \n")
printf (" or: rad = tb2rad2(temp_array, bandlist, inst='themis' etc)\n")
printf (" Input array of temps in partial set of 1 or more bands: temp_array[nx, ny, ibands] \n")
printf (" Broadband instruments only\n")
printf (" Dimension of ibands must match # of bands in bandlist \n")
printf (" Example: rad = tb2rad2(image, cat(3,5,7,x)) \n")
printf (" Returns rad[nx, ny, ibands] \n")
printf (" \n")
printf (" Note: Reads in appropriate radiance-to-temperature look-up table\n")
printf (" \n")
return(0)
}
# modified to include both 'tes=1' and 'inst="tes"' modes. prc `1/28/2011.
# changed from tb2rad2 to inst='lab1' mode, prc 12/16/2010
# if user sent in xaxis - assume they didn't also send in an 'inst'
if(HasValue(xaxis)) {
inst = 'spectrometer'
} else {
# check for 'old way' (tes=1) and convert to 'new way' (inst="tes")
args=cat("lab1", "lab2", "lab3", "tes", "tes5", "tes73", "mtes", "themis", "aster", "tims", "modis", axis=y)
for(i=1; i<=length(args); i+=1) {
if(HasValue(eval(args[,i]))==1) {
inst = args[,i]
}
}
}
# For spectrometers this routine accepts:
# 1 temp or an array of temps (z dim must equal # bands (e.g. 143 for TES))
# For broadband instrument this routine accepts:
# 1 temp, an array of temps, an array of temps in all bands for the instrument,
# or an array of temps in a subset of bands with a list of bands
# p. christenen 11/15/08 - derived from tbtorad
# tbtorad - p. christensen 10/2/01
# check for spectrometer or broadband
if(inst=='lab1' || inst=='lab2' || inst=='lab3' || inst=='tes' || inst=='tes5' || inst=='mtes' || inst=='spectrometer') {
# spectrometer
temp = $1
if(HasValue(xaxis)) {
x = xaxis
} else {
x = make_band(inst=inst)
}
nbands = dim(x)[3]
size=dim(temp)
rad = create(size[1], size[2], nbands, format=float)
# loop over x, y elements in input temp array
for(i=1; i<= size[1]; i+=1) {
for(j=1; j<= size[2]; j+=1) {
if(size[3] == 1) {
# sent in 1 temp
rad[i,j] = bbr(x, temp[i,j,1])
} else {
# sent in array of temps - 1 for each band
# z dimension must match nbands
if(size[3] != nbands) {
printf("z dimension (%d) of input temp array does not match nbands (%d)\n", size[3], nbands)
return(0)
} else {
rad[i,j] = bbr(x, temp[i,j])
}
}
}
}
# final test for radiances less than zero
rad[where(rad<0)] = 0.
return(rad)
} else {
# broadband instrument
# read in temp_to_radiance (temp_rad) array made by 'make_temp_rad' for broadband instruments
# temp_rad arrays are dimensioned: [xaxis+nbands, ntemps, 1]
# they contain the wavelength(microns) and nbands of integrated radiance values for each temperature
# default instrument is themis
# default path set up by Chris Edwards for davinci use on multiple platforms
PATH = $DV_SCRIPT_FILES
PATHI = PATH + "/instrument_parameters"
if(inst=='aster') {
temp_rad =read(PATHI+"/temp_rad_aster.vic")
} else if(inst=='tims') {
temp_rad =read(PATHI+"/temp_rad_tims.vic")
} else if(inst=='modis') {
# not done yet
} else if(inst=='themis') {
temp_rad = read(PATHI+'/temp_rad_v4')
} else {
# default is themis
temp_rad = read(PATHI+'/temp_rad_v4')
}
# set up the broadband instruments
nbands = dim(temp_rad)[1] - 1
ntemps = dim(temp_rad)[2]
if($ARGC == 1) {
# either 1 temp, an array of temps, or
# an nx x ny array in nbands, where nbands must match instrument
# no band list
temp = $1
size = dim(temp)
if(size[3] == 1) {
# 1 temp or an array of temps
if(size[1]== 1 && size[2] == 1) {
# correct size for 1 temperature
# clone temperature to match number of bands in temp_rad array
temp_array = clone(temp, 1, 1, nbands)
} else if(size[1] != 1 || size[2] != 1) {
# an array of temps with only 1 band and no list specified
temp_array = clone(temp,1, 1, nbands)
}
} else {
# sent in an array of temps that must be dimensioned [nx, ny, nbands]
if(size[3] != nbands) {
# a problem
printf("incorrect z-size - must match number in instrument\n")
return(0)
} else {
#printf("ok - %d by %d array with correct number of bands %d\n", size[1], size[2], size[3])
temp_array = temp
}
}
# using all the bands in the instrument - set temp_rad_interp to temp_rad
temp_rad_interp = temp_rad
} else if ($ARGC == 2) {
# sent in an array of temperatures in 1 or more bands and specified which band(s)
temp_array = $1
size = dim(temp_array)
bandlist = $2
list_size = dim(bandlist)
# check format of bandlist
if(list_size[2] >1 || list_size[3] > 1) {
printf("bad format for bandlist - bands should be in x; ny and nz should be 1\n")
return(0)
}
nlist = list_size[1]
if(size[3] != nlist) {
# a problem
printf("incorrect z-size - if specifing 1 or more bands, nbands in image must be same as nbands in list\n")
return(0)
}
if(nlist > nbands) {
# a problem
printf("incorrect band choice - higher than number of bands in this instrument\n")
return(0)
}
# test to make sure no band in the list is higher than # of bands
max_list = max(bandlist[,1,1])
if(max_list > nbands) {
# a problem
printf("a band in the list (%d) is higher than the number of bands (%d)\n", max_list, nbands)
return(0)
}
# ok - passed the tests
if(nlist == 1) {
# only using 1 band: move the correct band from temp_rad array into temp_rad_interp
temp_rad_interp = cat(temp_rad[1], temp_rad[1+bandlist], x)
} else {
# sent in a list of bands - move the appropriate bands of temp_rad into temp_rad_interp
temp_rad_interp = create(nlist+1, ntemps, 1, format=float)
temp_rad_interp[1] = temp_rad[1]
for(i=1; i<=nlist; i+=1) {
temp_rad_interp[i+1] = temp_rad[bandlist[i] + 1]
}
}
}
# now do the easy part - interpolate in radiance vs temperature array
nbands_used = dim(temp_rad_interp)[1]-1
printf("nbands=%d ntemps=%d size[1]=%d size[2]=%d size[3]=%d nbands_used=%d\n", nbands, ntemps, size[1], size[2], size[3], nbands_used)
rad = temp_array * 0.
for(i=1; i <= nbands_used; i+=1) {
index = i+1
rad[,,i] = interp(temp_rad_interp[index], temp_rad_interp[1], temp_array[,,i])
}
# clean up to free memory
temp_array = 0.
# final test for radiances less than zero
rad[where(rad<0)] = 0.
return(rad)
}
}
define make_band(inst, cm, wave, full, descend, apod, lab1, lab2, lab3, tes, tes5, tes73, mtes, themis, aster, tims, modis, master) {
# makes wavelength centers (and min and max) in wavelength and wavenumber
# for IR instruments
# instruments supported are lab1, lab2, lab3, tes, tes5, tes73, themis, aster, tims, master
# help section is at the bottom of the code block
# added master to instruments supported. 4/29/11 p. christensen
# modified to include both 'tes=1' and 'inst="tes"' modes. Name changed back to make_band() prc `1/28/2011.
# changed to 'inst= "_"' mode 1/22/2011 prc. name changed to make_band3()
# p. christensen 11/28/08
# default path set up by Chris Edwards for davinci use on multiple platforms
PATH = $DV_SCRIPT_FILES
PATHI = PATH +"/instrument_parameters"
struc = struct()
# do initial tests
# check for 'old way' (tes=1) and convert to 'new way' (inst="tes)
args=cat("lab1", "lab2", "lab3", "tes", "tes5", "tes73", "mtes", "themis", "aster", "tims", "modis", "master", axis=y)
for(i=1; i<=length(args); i+=1) {
if(HasValue(eval(args[,i]))==1) {
inst = args[,i]
}
}
if(HasValue(wave) && HasValue(cm)) {
printf("Can't request both 'wave' and 'cm'\n")
return(0)
}
if(HasValue(full) && HasValue(descend)) {
printf("Please don't request 'descend' in 'full' mode. Keep life simple\n")
return(0)
}
# check the instruments
if(inst == "lab1" || inst == "lab2" || inst == "lab3" || inst == "tes" || inst == "tes5" || inst == "tes73" || inst == "mtes") {
# do the spectrometers
if(HasValue(full) && HasValue(wave)) {
printf("'full' mode only returns default mode - don't use 'wave' for spectrometers or 'descend'\n")
return(0)
}
if(HasValue(cm) == 0 && HasValue(wave) == 0) {
# default for spectrometers is wavenumber
cm = 1
}
if(inst == "tes" || inst == "tes73") {
# borrowed from make_tesx prc 2/3/09
npts = 143
first_pt = 14
# set prime factor to default for single scan (1344; Det. 2 & 5)
nprime = 1344
neon_wavelength = 0.7032e-4
tes_res = 1. / (nprime * neon_wavelength)
# old make_xt2 hardcoded value tes_res = 10.5808820
bandcenter_cm = create(1,1,npts, start = tes_res * first_pt, step = tes_res, format=float)
# bandcenter_cm = make_tesx()
# truncate after 4 decimal places - legacy
bandcenter_cm = bandcenter_cm * 10000.
bandcenter_cm = int(bandcenter_cm)
bandcenter_cm = bandcenter_cm/10000.
nwave = dim(bandcenter_cm)[3]
if(HasValue(apod)) {
# do the approximation to S. Chase's apodization for center detector 2
# Data are from TES II Calibration Report Appendix 11/12/98
# subtract 0.5 to make the slope exactly half that of 5cm-1 case and give
# consistent results with 5cm-1 case and Calib Report Appendix
slope = (5.27 - 0.438) / (nwave-.5)
offset = .438
for(k=1; k<=nwave; k+=1) {
bandcenter_cm[1,1,k] = bandcenter_cm[1,1,k] + ((k-1)* slope + offset)
}
}
#make a 73 point tes spectrum
if(inst == "tes73") {
bandcenter_cm=cat(bandcenter_cm[,,9:35],bandcenter_cm[,,65:110],axis=z)
}
} else if(inst == "tes5") {
npts = 143 * 2
first_pt = 28
# set prime factor to default for double scan (1344*2; Det. 2 & 5)
nprime = 1344 * 2
# tes resolution for 1344 pt fft
neon_wavelength = 0.7032e-4
tes_res = 1. / (nprime * neon_wavelength)
# old make_xt2 hardcoded value tes_res = 10.5808820 / 2.
bandcenter_cm = create(1,1,npts, start = tes_res * first_pt, step = tes_res, format=float)
# bandcenter_cm = make_tesx(scan_len=2)
# truncate after 4 decimal places - legacy
bandcenter_cm = bandcenter_cm * 10000.
bandcenter_cm = int(bandcenter_cm)
bandcenter_cm = bandcenter_cm/10000.
nwave = dim(bandcenter_cm)[3]
if(HasValue(apod)) {
# do the approximation to S. Chase's apodization for center detector 2
# Data are from TES II Calibration Report Appendix 11/12/98
slope = (5.27 - 0.438) / (nwave-1)
offset = .438
for(k=1; k<=nwave; k+=1) {
bandcenter_cm[1,1,k] = bandcenter_cm[1,1,k] + ((k-1)* slope + offset)
}
}
} else if(inst == "mtes") {
# copied from make_mtesx()
npts = 512
start_pt = 35
end_pt = 201
# set wavenumber of first point downlinked (14 or 28 in PROM) #
# set prime factor
nprime = 1024;
laser_wavelength = 0.978
delta_x = (1./(nprime * laser_wavelength))*10000.
mtx = create(npts, 1, 1, format=float)
for(i=1; i<= npts; i+=1) {
# index (i) of 1 = DC term in FFT; index 2 = 9.9853 cm-1; etc
# index 1 = channel 0; etc
mtx[i] = float(i-1) * delta_x
}
# truncate after 4 decimal places
mtx = mtx * 10000.
mtx = int(mtx)
mtx = mtx/10000.
mtx = translate(mtx,x,z)[,,start_pt:end_pt]
bandcenter_cm = mtx
# bandcenter_cm = make_mtesx()
} else if(inst == "lab1") {
# compute the values
# taken from best fit to xaxis spectrum from Nicholet
first_pt = 221.777
# derived from best fit to xaxis spectrum from Nicholet
lab_res = 1.9284983
npts = 923
bandcenter_cm = create(1,1,npts, start = first_pt, step = lab_res, format=float)
# read in the values from a spectrum from the Nicholet
# bandcenter_cm = read(PATHI+"/lab1cmx.vic")
} else if(inst == "lab2") {
# compute the values
# taken from best fit to xaxis spectrum from Nicholet
first_pt = 200.564
# derived from best fit to xaxis spectrum from Nicholet
lab_res = 1.9284983
npts = 1971
bandcenter_cm = create(1,1,npts, start = first_pt, step = lab_res, format=float)
# read in the values from a spectrum from the Nicholet
# bandcenter_cm = read(PATHI+"/lab2cmx.vic")
} else if(inst == "lab3") {
# compute the values
# taken from best fit to xaxis spectrum from Nicholet
first_pt = 380.001
# derived from best fit to xaxis spectrum from Nicholet
lab_res = 1.9284983
npts = 840
bandcenter_cm = create(1,1,npts, start = first_pt, step = lab_res, format=float)
# read in the values from a spectrum from the Nicholet
# bandcenter_cm = read(PATHI+"/lab2cmx.vic")
}
# convert to wavelength in case it is wanted
bandcenter = cm2micron(bandcenter_cm)
# return all values if requested - only do default modes for 'full'
if(HasValue(full)) {
struc.bandcenter = bandcenter
struc.bandcenter_cm = bandcenter_cm
return(struc)
}
if(HasValue(wave)) {
# return wavelength
if(HasValue(descend)) {
# when doing spectrometers, wavelength version is already descending
} else {
# make wavelength into ascending order (default)
bandcenter = reverse_spectrum(bandcenter)
}
return(bandcenter)
} else {
# return wavenumber (default)
if(HasValue(descend)) {
# make wavenumber into descending order
bandcenter_cm = reverse_spectrum(bandcenter_cm)
} else {
# wavelength version is already ascending (default)
}
return(bandcenter_cm)
}
} else if(inst == "themis" || inst == "aster" || inst == "tims" || inst == "modis" || inst == "master") {
# do the broadband instruments
if(HasValue(full) && HasValue(cm)) {
printf("'full' mode only returns default mode - don't use 'cm' for broadband inst. or 'descend'\n")
return(0)
}
if(HasValue(cm) == 0 && HasValue(wave) == 0) {
# default for broadband instruments is wavelength
wave = 1
}
if (inst == "themis") {
# Make 10-point THEMIS wavelength and wavenumber axes
# source THEMIS calibration report
bandcenter=create(1,1,10, format=float)
bandmin=create(1,1,10, format=float)
bandmax=create(1,1,10, format=float)
bandcenter[1,1,1] = 6.78
bandcenter[1,1,2] = 6.7801
bandcenter[1,1,3] = 7.93
bandcenter[1,1,4] = 8.56
bandcenter[1,1,5] = 9.35
bandcenter[1,1,6] = 10.21
bandcenter[1,1,7] = 11.04
bandcenter[1,1,8] = 11.79
bandcenter[1,1,9] = 12.57
bandcenter[1,1,10] = 14.88
# source: P. Christensen analysis of themis_IR_filter_response_norm_v2.davinci
# wavelength is half-power (50%) point of normalized response
bandmin[1,1,1] = 6.277
bandmin[1,1,2] = 6.2775
bandmin[1,1,3] = 7.382
bandmin[1,1,4] = 7.984
bandmin[1,1,5] = 8.756
bandmin[1,1,6] = 9.668
bandmin[1,1,7] = 10.455
bandmin[1,1,8] = 11.265
bandmin[1,1,9] = 12.175
bandmin[1,1,10] = 14.46
# source: P. Christensen analysis of themis_IR_filter_response_norm_v2.davinci
# wavelength is half-power (50%) point of normalized response
bandmax[1,1,1] = 7.285
bandmax[1,1,2] = 7.2855
bandmax[1,1,3] = 8.470
bandmax[1,1,4] = 9.140
bandmax[1,1,5] = 9.950
bandmax[1,1,6] = 10.755
bandmax[1,1,7] = 11.635
bandmax[1,1,8] = 12.325
bandmax[1,1,9] = 12.987
bandmax[1,1,10] = 15.32
} else if(inst == "aster") {
# Make 5-point ASTER wavelength and wavenumber axes
# source: www.science.aster.ersdac.or.jp/en
bandcenter=create(1,1,5, format=float)
bandmin=create(1,1,5, format=float)
bandmax=create(1,1,5, format=float)
bandcenter[1,1,1] = 8.274
bandcenter[1,1,2] = 8.626
bandcenter[1,1,3] = 9.072
bandcenter[1,1,4] = 10.654
bandcenter[1,1,5] = 11.303
# source: P. Christensen analysis of aster_ir_spectral_response_wave.vic
# wavelength is half-power (50%) point of normalized response
bandmin[1,1,1] = 8.105
bandmin[1,1,2] = 8.455
bandmin[1,1,3] = 8.895
bandmin[1,1,4] = 10.325
bandmin[1,1,5] = 11.005
# source: P. Christensen analysis of aster_ir_spectral_response_wave.vic
# wavelength is half-power (50%) point of normalized response
bandmax[1,1,1] = 8.45
bandmax[1,1,2] = 8.800
bandmax[1,1,3] = 9.255
bandmax[1,1,4] = 10.985
bandmax[1,1,5] = 11.595
} else if (inst == "tims") {
# Make 6 point TIMS wavelength and wavenumber axes
# source: PRC analyis of tims.resp files from Elsa Abbot 12/18/08
# Made by 'make_tims_resp' - average of 10 tims.resp files
bandcenter=create(1,1,6, format=float)
bandmin=create(1,1,6, format=float)
bandmax=create(1,1,6, format=float)
bandcenter[1,1,1] = 8.38450
bandcenter[1,1,2] = 8.78470
bandcenter[1,1,3] = 9.18000
bandcenter[1,1,4] = 9.89410
bandcenter[1,1,5] = 10.7344
bandcenter[1,1,6] = 11.6222
# source: P. Christensen analysis of tims_ir_spectral_response_avg_wave.vic
bandmin[1,1,1] = 8.17
bandmin[1,1,2] = 8.58
bandmin[1,1,3] = 9.00
bandmin[1,1,4] = 9.53
bandmin[1,1,5] = 10.29
bandmin[1,1,6] = 11.28
bandmax[1,1,1] = 8.55
bandmax[1,1,2] = 8.95
bandmax[1,1,3] = 9.31
bandmax[1,1,4] = 10.21
bandmax[1,1,5] = 11.18
bandmax[1,1,6] = 11.79
} else if (inst == "modis") {
} else if (inst == "master") {
# thermal ir bands only
# Make 5-point MASTER wavelength and wavenumber axes
# source: c. edwards 4/28/11
bandcenter=create(1,1,10, format=float)
bandmin=create(1,1,10, format=float)
bandmax=create(1,1,10, format=float)
# source: P. Christensen analysis (by finding max of normaized bandpasses)
# of master_ir_spectral_response_avg_wave.vic
bandcenter[1,1,1] = 7.793
bandcenter[1,1,2] = 8.215
bandcenter[1,1,3] = 8.631
bandcenter[1,1,4] = 9.093
bandcenter[1,1,5] = 9.716
bandcenter[1,1,6] = 10.131
bandcenter[1,1,7] = 10.647
bandcenter[1,1,8] = 11.311
bandcenter[1,1,9] = 12.177
bandcenter[1,1,10] = 12.918
# source: P. Christensen hand analysis (by plotting) of master_ir_spectral_response_avg_wave.vic
bandmin[1,1,1] = 7.59
bandmin[1,1,2] = 7.93
bandmin[1,1,3] = 8.44
bandmin[1,1,4] = 8.89
bandmin[1,1,5] = 9.50
bandmin[1,1,6] = 9.93
bandmin[1,1,7] = 10.38
bandmin[1,1,8] = 11.05
bandmin[1,1,9] = 11.91
bandmin[1,1,10] = 12.68
bandmax[1,1,1] = 7.92
bandmax[1,1,2] = 8.40
bandmax[1,1,3] = 8.80
bandmax[1,1,4] = 9.29
bandmax[1,1,5] = 9.92
bandmax[1,1,6] = 10.34
bandmax[1,1,7] = 10.99
bandmax[1,1,8] = 11.72
bandmax[1,1,9] = 12.45
bandmax[1,1,10] = 13.15
} else {
printf("No instrument specified\n")
return(0)
}
# convert to wavenumber in case it is wanted
bandcenter_cm = 10000./bandcenter
bandmin_cm = 10000./bandmin
bandmax_cm = 10000./bandmax
# return all values if requested - only do default modes for 'full'
if(HasValue(full)) {
struc.bandcenter = bandcenter
struc.bandcenter_cm = bandcenter_cm
struc.bandmin = bandmin
struc.bandmin_cm = bandmin_cm
struc.bandmax = bandmax
struc.bandmax_cm = bandmax_cm
return(struc)
}
if(HasValue(cm)) {
if(HasValue(descend)) {
# when doing braodband, wavenumber version is already descending
} else {
# make wavenumber into ascending order (default)
bandcenter_cm = reverse_spectrum(bandcenter_cm)
}
return(bandcenter_cm)
} else {
# return wavelength (default)
if(HasValue(descend)) {
# make wavelength into descending order
bandcenter = reverse_spectrum(bandcenter)
} else {
# wavenumber version is already ascending (default)
}
return(bandcenter)
}
} else {
printf("No valid instrument specified\n\n")
printf (" Make bands ('xaxes') for major IR instruments\n")
printf (" Usage: x = make_band([inst='inst'] or [tes=1][tes5=1][tes73=1][mtes=1][themis=1][lab1=1][lab2=1]lab3=1][aster=1][tims=1][modis=1][full=1][cm=1][wave=1][descend=1]) \n")
printf (" Default is to return wavelength or wavenumber in natural order for specified instrument\n")
printf (" works in inst='inst' (e.g. inst='themis') or themis=1 modes\n")
printf (" i.e. ascending wavenumbers for spectrometers; ascending wavelength for broadband instruments\n")
printf (" [full = 1] Returns bandcenter and min and max (where appropriate) in cm-1 and microns\n")
printf (" 'full' mode only returns default mode (no descend)\n")
printf (" [cm = 1] Returns bandcenter in ascending cm-1\n")
printf (" [wave = 1] Returns bandcenter in ascending wavelength\n")
printf (" [descend = 1] Inverts default order\n")
printf (" [apod = 1] Applies approximation of S. Chase apodization to TES\n")
printf (" If cm flag is specified for broadband or wave flag for spectrometers, the order is still ascending\n")
printf (" descend flag inverts the default order\n")
printf (" \n")
printf (" p. christensen 11/28/08\n")
return(0)
}
}
define make_temp_rad(themis, aster, tims) {
if (HasValue(themis)==0 && HasValue(aster)==0 && HasValue(tims)==0) {
printf (" make tbtorad and radtotb look up table by integrating Planck function times THEMIS IR filter functions\n")
printf (" \n")
printf (" default is THEMIS, also does ASTER and TIMS \n")
printf (" \n")
printf (" creates array of temperature and spectral radiance, weighted by filter response \n")
printf (" spectral radiance is in units of W cm-2 str-1 micron-1 \n")
printf (" Array is [nbands+1,ntemps], where x dimension is temp and n spectral radiances \n")
printf (" [themis = 1] \n")
printf (" [aster = 1] \n")
printf (" [tims = 1] \n")
return(0)
}
# p.christensen 10/2/01
# generalized to include ASTER, TIMS, and MODIS: prc 11/08
######## history
# THEMIS is still version 4
# P. Christensen 11/14/08
# added PATH
# changed to spectral radiance - p. christensen 3/29/02
# increment to version 4
# changed start and end temps to 100. and 380. - p. christensen 3/14/02
# increment to version 3
#start_temp = 120.
#end_temp = 360.
# used up to 11/10/01 - and beyond
#d_temp = 0.1
# experiment - prc - 11/11/01
# only changed final temps by <0.001 prc
# wrote results to /themis/data/archive/temp_rad_5
# tested in test_temp_rad5
# used d_temp = 0.5 from here on
# v2 - prc 1/4/02
# used better relative spectral response - filter * window * beamsplitter (/themis/calib/themis_IR_response_norm_v2.davinci)
# used d_temp = 0.5
#################
d_temp = 0.5
start_temp = 100.
end_temp = 380.
# set path
# default path set up by Chris Edwards for davinci use on multiple platforms
PATH = $DV_SCRIPT_FILES
PATHI = PATH +"/instrument_parameters"
ntemps = int( (end_temp - start_temp) / d_temp + 1 )
# read in filter response functions [nbands+1 x nwaves x 1] - x[1] = wavelengths; x[2-n] = responses for each filter
# for themis f=[11,4000,1] - 1 xaxis + 10 bands at 4000 wavelengths
if(HasValue(aster)) {
f = read(PATHI+"/aster_ir_spectral_response_wave.vic")
} else if(HasValue(tims)) {
f = read(PATHI+"/tims_ir_spectral_response_avg_wave.vic")
} else if(HasValue(modis)) {
# not done yet
} else if(HasValue(themis)) {
f = ascii(PATHI+"/themis_IR_filter_response_norm_v2.davinci",format=float)
} else {
f = ascii(PATHI+"/themis_IR_filter_response_norm_v2.davinci",format=float)
}
nbands = int(dim(f)[1] - 1)
nwaves = int(dim(f)[2])
printf("start temp = %.1f end temp = %.1f d temp = %.1f\n", start_temp, end_temp, d_temp)
printf("nwaves = %d ntemps = %d nbands = %d\n", nwaves, ntemps, nbands)
# create array for spectral radiance for each band, each temperature and set to 0.
rsum = create(nbands, ntemps, 1, format=double)
rsum = rsum * 0.
# create array for temperatures
temps = create(1,ntemps,1,format=double)
temps = temps * 0.
# make array of delta wave
dwave = create(1, nwaves, format=double)
for(i=1; i<= nwaves-1; i+=1) {
dwave[,i] = f[1,i+1]-f[1,i]
}
dwave[,nwaves] = dwave[,nwaves-1]
# make bandpass integral
# sum filter response over all wavelengths
fsum = sum((f[2:(nbands+1)] * dwave),y)
# loop on temps
for(i=1; i<= ntemps; i+=1) {
# compute temp for this index
temp = start_temp + (i-1) * d_temp
# load temp array
temps[1,i,1]=temp
# get spectral radiance (W cm-2 micron-1 str-1) for each wavelength for this temperature
# f[1,,] contains wavelengths
rad = bbrw(f[1,,], temp)
# multiply spectral radiance * filter_response * dwave (micron) and sum over all wavelengths
# result is spectral radiance - W cm-2 str-1 micron-1
rsum[,i,] = sum( (rad * f[2:(nbands+1)] * dwave) ,y)
}
# compute spectral radiance (units of W cm-2 str-1 micron-1) - added prc 3/29/02
sp_rsum = rsum/fsum
# create final temp_rad array
temp_rad = temps//sp_rsum
#temp_rad = temps//float(sp_rsum)
if (HasValue(write)) {
# verson 4 3/12902 - spectral radiance
if(HasValue(aster)) {
name = "/temp_rad_aster.vic"
fileout = PATHI+ name
printf("writing file to %s\n", fileout)
write(temp_rad, fileout, vicar, force=1)
} else if(HasValue(tims)) {
fileout = PATHI+"/temp_rad_tims.vic"
write(temp_rad, fileout, vicar, force=1)
} else if (HasValue(modis)) {
# not done yet
} else {
# themis is default
# verson 4 3/12902 - spectral radiance
fileout = PATHI+"/temp_rad_v4"
printf("Writing is disabled for THEMIS to avoid accidentally changing\n")
printf(" the official THEMIS calibration files - PRC 12/19/08\n")
#printf("Writing temp_out file %s\n", fileout)
#write(temp_rad, fileout, vicar, force=1)
#write(fsum,"/themis/calib/ir_band_width_v2",vicar,force=1)
#write(fsum,"/themis/data/archive/ir_band_width_v2",vicar,force=1)
}
# verson 2 1/4/02 - new filter * window * beamsplitter data
# fileout = "/themis/calib/temp_rad_v2"
# printf("Writing temp_out file %s\n", fileout)
# write(temp_rad, fileout, vicar, force=1)
# write(temp_rad,"/themis/data/archive/temp_rad_v2",vicar,force=1)
# write(temp_rad,"/themis/data/archive/temp_rad_v1.1",vicar,force=1)
# write(temp_rad,"/themis/data/archive/temp_rad",vicar,force=1)
# write(fsum,"/themis/calib/ir_band_width_v2",vicar,force=1)
# write(fsum,"/themis/data/archive/ir_band_width_v2",vicar,force=1)
}
return(temp_rad)
}
define i2i(to, from, wave, cm, full, filters, descend, list, supported) {
if ($ARGC == 0 && HasValue(supported)==0 ) {
printf(" Convert spectra from one IR instrument to another\n")
printf(" \n")
printf(' Usage: spec_out = i2i(spec, from="from_inst", to="to_inst" [wave=1][cm][full][filters][descend])\n')
printf(" Allowed 'from' instruments are lab1, lab2, tes, tes5, tes73, and mtes\n")
printf(" Allowed 'to' instrument are lab1, lab2, tes, tes5, tes73, mtes, themis, aster, tims, master\n")
printf(" Both 'from' and 'to' can be arrays of xaxis values. Size of 'from' must equal size of spec\n")
printf(" Use 'to'=xaxis or 'from'=xaxis' to create or use customized output or input spectra.\n")
printf(" Default instrument input (from) spectra must be in ascending, cm-1 form\n")
printf(" Exception is converting broadband instrument from ascending wavelength\n")
printf(" to ascending cm\n")
printf(" This mode also returns the xaxis to avoid any confusion\n")
printf(" Ex: a = i2i(spect, from='themis', to='themis') returns 'a' in ascending cm form\n")
printf(" Default 'to' spectra are in ascending (low to high) cm for spectrometers\n")
printf(" or ascending wavelength for broadband instruments\n")
printf(" [wave=1] returns spectra in ascending wavelength form\n")
printf(" [cm=1] returns spectra in ascending wavenumber (cm) form\n")
printf(" [descend=1] returns spectrum in descending wavelength or wavenumber order\n")
printf(" Note: Non-standard cm or wave modes return both spectra and xaxis\n")
printf(" [list = list]. In spectral library mode, only move a portion of the resampled data\n")
printf(" Ex: list = '1//4//7//21:27, 33//35//51:60'\n")
printf(" [full=1] returns both spectrum and xaxis\n")
printf(" [filters=1] returns spectrum, xaxis, and filter functions\n")
printf(" \n")
printf(" Examples: speca = i2i(spec, from = 'lab1', to ='themis')\n")
printf(" speca = i2i(spec, from = 'lab1', to ='aster', wave=1)\n")
printf(" speca = i2i(spec, from = 'lab1', to ='lab1', wave=1)\n")
printf(" speca = i2i(spec, from = 'tes', to = 'aster', cm=1)\n")
printf(" speca = i2i(spec, from = 'lab1', to = my_out_axis)\n")
printf(" speca = i2i(spec, from = my_in_axis, to = my_out_axis)\n")
printf(" \n")
printf(" Spectral library mode examples:\n")
printf(" speclib = i2i(asulib, to ='aster', list='1//3//21:28, //33//41')\n")
printf(" \n")
printf(" p.christensen 11/08, 7/10\n")
return(0)
}
# added master to list of allow 'to' instrument
# modified to resample spectral library structures and move all elements,
# including sublib option using 'list' - prc 4/25/11
# modified to do special cases of from = to - prc
# reorganized significantly to not use to and from internally - prc 7/1/10
# changed to use 'to' and 'from' keywords - prc 6/30/2010
# added custom and xtospec - prc 3/24/2010
# removed resample option from broadband instruments - prc 1/18/09
# added list of supported insturments and TES 73 point spectra
#HERE is the description list of supported instruments
if(HasValue(supported)) {
list={}
list.list=cat("lab1",\
"lab2",\
"tes",\
"tes5",\
"tes73",\
"mtes",\
"themis",\
"tims",\
"aster",\
"master",axis=y);
list.description=cat("SpecLab1 sampling (200-2000 @ 2cm-1)",\
"SpecLab2 sampling (200-4000 @ 2cm-1)",\
"TES single sampling (148-1650 @ 10cm-1)",\
"TES double sampling (148-1655 @ 5cm-1)",\
"TES 73-point sampling (339-1301 @ 10cm-1)",\
"Mini-TES sampling (339-1997 @ 10cm-1)",\
"THEMIS IR sampling",\
"TIMS IR sampling",\
"ASTER IR sampling",\
"MASTER IR sampling",axis=y);
return(list);
}
# set up
in = $1
# test to see if user sent in a single standard spectral structure
if (type(in) == "STRUCT") {
if(HasValue(get_struct(in, "data")) == 0) {
printf("You sent in a structure, but it doesn't have a 'data' field\n")
return(0)
}
if(HasValue(from)==0){
from = in.xaxis
}
# set structure logical flag
st = 1
spec = in.data
} else {
# user sent in a spectrum
st = 0
spec = in
}
struc = struct()
# set path
# default path set up by Chris Edwards for davinci use on multiple platforms
PATH = $DV_SCRIPT_FILES
PATHI = PATH + "/instrument_parameters"
if(HasValue(filters)) {
# automaticaly set mode to 'full' if user requested the filter functions
full = 1
}
# do the initial checking of to and from, keywords or arrays
if(HasValue(from) == 0) {
printf("No from instrument specified\n")
return(0)
}
if(HasValue(to) == 0) {
printf("No to instrument specified\n")
return(0)
}
# check to see if 'to' is an x axis
format = format(to)
if(format != 'STRING') {
size = dim(to)
xto = to
inst_to = 'custom_to'
} else {
inst_to = to
}
# check to see if 'from' is an x axis
format = format(from)
if(format != 'STRING') {
size = dim(from)
xfrom = from
inst_from = 'custom_from'
} else {
inst_from = from
}
# check for special cases; these are ok only if doing the special case of
# converting broadband instrument from ascending wavelength to ascending wavenumber
if(inst_from == "themis" || inst_from == "aster" || inst_from == "tims" || inst_from == "modis"|| inst_from == "master") {
if(inst_to == "themis" && inst_from == "themis") {
xto = make_band(themis=1, cm=1)
} else if(inst_to == "aster" && inst_from == "aster") {
xto = make_band(aster=1, cm=1)
} else if(inst_to == "tims" && inst_from == "tims") {
xto = make_band(tims=1, cm=1)
} else if(inst_to == "master" && inst_from == "master") {
xto = make_band(master=1, cm=1)
} else {
printf("themis, aster, tims, modis, and master can only 'inst_from' if also 'to'\n")
return(0)
}
# invert spectrum
data = reverse_spectrum(spec)
if(HasValue(full)) {
# this modes does not return filter functions
struc.data = data
struc.xto = xto
return(struc)
} else {
return(data)
}
}
# check for more special cases; these are for doing the special case of
# converting spectrometer instrument from ascending wavenumber to ascending wavelength
if(inst_to == inst_from) {
if(inst_to == 'lab1' && inst_from == 'lab1') {
xto = make_band(lab1=1, wave=1)
} else if(inst_to == 'lab2' && inst_from == 'lab2') {
xto = make_band(lab2=1, wave=1)
} else if(inst_to == 'tes' && inst_from == 'tes') {
xto = make_band(tes=1, wave=1)
} else if(inst_to == 'tes5' && inst_from == 'tes5') {
xto = make_band(tes5=1, wave=1)
} else if(inst_to == 'tes73' && inst_from == 'tes73') {
xto = make_band(tes73=1, wave=1)
} else if(inst_to == 'mtes' && inst_from == 'mtes') {
xto = make_band(mtes=1, wave=1)
}
# invert spectrum
data = reverse_spectrum(spec)
if(HasValue(full)) {
# this modes does not return filter functions
struc.data = data
struc.xto = xto
return(struc)
} else {
return(data)
}
}
# do remaining checks
# check for legitimate 'from' instruments' for 'normal' cases
if(inst_from == "tes" || inst_from == "tes5" || inst_from == "tes73" || inst_from == "mtes" || inst_from == "lab1" || inst_from == "lab2" || inst_from == "custom_from") {
# these are ok - go on
} else {
printf("incorrect input instrument\n")
printf("Only lab1, lab2, tes, tes5, and mtes allowed for 'from' at this time\n")
printf("except in special broadband case of 'to' = 'from'\n")
return(0)
}
# check for legitimate 'to' instruments'
if(inst_to == "lab1" || inst_to == "lab2" || inst_to == "tes" || inst_to == "tes5" || inst_to == "mtes" || inst_to == "themis" || inst_to == "aster" || inst_to == "tims" || inst_to == "modis" || inst_to == "master" || inst_to=="tes73" || inst_to == "custom_to") {
} else {
printf("incorrect output instrument\n")
return(0)
}
# get the xaxis for 'from' spectra
if(HasValue(xfrom)==0) {
# no xaxis sent in (default) - get the 'from' xaxis
# xaxes are in cm form, ascending (low cm to high cm)
if(inst_from == "lab1") {
xfrom = make_band(lab1=1)
} else if(inst_from == "lab2") {
xfrom = make_band(lab2=1)
} else if(inst_from == "tes") {
xfrom = make_band(tes=1)
} else if(inst_from == "tes5") {
xfrom = make_band(tes5=1)
} else if(inst_from == "tes73") {
xfrom = make_band(tes73=1)
} else if(inst_from == "mtes") {
xfrom = make_band(mtes=1)
} else {
printf("Only lab1, lab2, tes, tes5, tes73, and mtes allowed for 'from' at this time\n")
return(0)
}
}
size = dim(spec)
nxfrom = size[1]
nyfrom = size[2]
# spec_len is length of from spectrum
spec_len = size[3]
# xfrom_len is length of 'from' 'xaxis'
xfrom_len = dim(xfrom)[3]
# check to make sure length of from spectra and xaxis match
if(xfrom_len != spec_len) {
printf("Length of spectrum (%d) does not match that of xaxis (%d)\n", spec_len, xfrom_len)
return(0)
}
# work through possible 'to' combinations
if(inst_to == "lab1" || inst_to == "lab2" || inst_to == "tes" || inst_to == "tes5" || inst_to == "mtes" || inst_to=="tes73" || inst_to == "custom_to") {
# the spectrometers - only do the resample approach
# xto is in cm form, ascending
if(inst_to == "lab1") {
xto = make_band(lab1=1)
} else if(inst_to == "lab2") {
xto = make_band(lab2=1)
} else if(inst_to == "tes") {
xto = make_band(tes=1)
} else if(inst_to == "tes5") {
xto = make_band(tes5=1)
} else if(inst_to == "tes73") {
xto = make_band(tes73=1)
} else if(inst_to == "mtes") {
xto = make_band(mtes=1)
} else if(inst_to == "custom_to") {
# do nothing - xaxis already set to xto
}
# do the resampling - ok for hyperspectral spectrometers
if($DV_OS=="win") {
specnew = resampleudf(spec, xfrom, xto).data
} else {
specnew = resample(spec, xfrom, xto).data
}
} else {
printf("\ndoing broadband\n")
# the broadband instruments - do a convolue
# read in spectral response versus wavelength
# wavenumber and response function are in ascending wavenmuber order
if(inst_to == "themis") {
# make xto in cm form, ascending (low cm to high cm) (high wave to low wave)
xto = make_band(themis=1, cm=1)
fati=read(PATHI+"/themis_ir_spectral_response_cm.vic")
} else if(inst_to == "aster") {
# make xto in cm form, ascending (low cm to high cm) (high wave to low wave)
xto = make_band(aster=1, cm=1)
fati=read(PATHI+"/aster_ir_spectral_response_cm.vic")
} else if(inst_to == "tims") {
# make xto in cm form, ascending (low cm to high cm) (high wave to low wave)
xto = make_band(tims=1, cm=1)
fati=read(PATHI+"/tims_ir_spectral_response_avg_cm.vic")
} else if(inst_to == "master") {
xto = make_band(master=1, cm=1)
fati=read(PATHI+"/master_ir_spectral_response_avg_cm.vic")
} else if(inst_to == "modis") {
}
# convolue input spectrum with instrument spectral response
# make sure input (from) xaxis is in ascending order
if(xfrom[1,1,xfrom_len] <= xfrom[1,1,1]) {
printf("input xaxis is in descending wavenumber range\n")
return(0)
}
# first resample spectral response array to input spectrum
# nx is number bands +1(xaxis) in spectral response array
nx= dim(fati)[1]
nbands = nx - 1
# nwaves is length of 'to' 'xaxis'
nwaves = dim(xto)[3]
# create output array
specnew = create(nxfrom, nyfrom, nbands, format=float)
# create new (resampled) spectral response array
filt_resample = create(nbands, 1, xfrom_len, format=float)
# loop over bands
for(i=1; i<=nbands; i+=1) {
# resample response function of 'to' instrument to 'from' xaxis
if($DV_OS=="win") {
resampled_response = resampleudf(fati[i+1], fati[1], xfrom)
} else {
resampled_response = resample(fati[i+1], fati[1], xfrom)
}
filt_resample[i] = resampled_response.data
specnew[,,i] = sum(spec * filt_resample[i], axis=z) / sum(filt_resample[i],axis=z)
}
# input spectrum was in ascending cm order; spectral response is in
# ascending cm order, and xto is in ascending cm order,
# but specnew is in band order - i.e. ascending
# wavelength order, descending cm order. Must invert specnew to put it into
# ascending cm order
specnew = reverse_spectrum(specnew)
if(HasValue(filters)) {
printf("set up to return filter functions\n")
struc.filters = filt_resample
struc.xfilters = xfrom
}
# set default condition for returned spectrum to 'wavelength' unless user specified 'cm'
if(HasValue(cm) == 0){
wave = 1
}
}
if(HasValue(wave)) {
# convert 'to' spectra to ascending wavelength form'
xto = cm2micron(xto)
# invert xaxis
xto = reverse_spectrum(xto)
# invert spectrum
specnew = reverse_spectrum(specnew)
if(HasValue(filters)) {
# invert filter function
struc.filters = reverse_spectrum(struc.filters)
# convert xfilters to wavelength and invert
struc.xfilters = reverse_spectrum(cm2micron(struc.xfilters))
}
}
if(HasValue(descend)) {
# invert both spectra and xaxis
specnew = reverse_spectrum(specnew)
xto = reverse_spectrum(xto)
}
# return data
if(st == 1) {
# send back a structure with new data and xaxis
if(HasValue(list)==0) {
# move entire structure
out_struc = in
# replace data and xaxis with resampled values
out_struc.data = specnew
out_struc.xaxis = xto
return(out_struc)
} else {
# move in user-specified vectors from structure
in.data = specnew
out_struc = sublib(in, list)
# reset the xaxis with the resampled xaxis
out_struc.xaxis = xto
return(out_struc)
}
} else {
if(HasValue(full)) {
struc.data = specnew
struc.xaxis = xto
return(struc)
} else {
# only return spectrum (default)
return(specnew)
}
}
}
define i2i_tutorial(index) {
printf("read in ASU spectral library\n")
cmd = "asu = read($DV_EX+'/ASU_minlib.hdf')"
echo(cmd)
eval(cmd)
pause("pausing")
if(HasValue(index)) {
i = index
} else {
i = 1
}
printf("\n**** Resample a lab spectrum to ASTER resolution and plot both vs wavenumber **\n\n")
printf("create lab1 x axis (default is wavenumber mode for spectrometers)\n")
cmd = "xlab1 = make_band(inst='lab1')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("create aster x axis in wavenumber mode\n")
cmd = "xastercm = make_band(inst='aster', cm=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample 1 lab spectrum to ASTER resolution with wavenumber xaxis\n")
cmd = "aster = i2i(asu.data[i], from = 'lab1', to = 'aster', cm=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("plot lab and resampled ASTER spectra\n")
cmd = "pplot({asu.data[i], aster}, {'Lab Spectrum', 'ASTER'}, {1,5}, xaxis={xlab1, xastercm}, xlabel = 'Wavenumber', ylabel = 'Emissivity', key = '1400, .95')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("\n*** A simpler approach ***\n\n")
printf("Resample to ASTER and return xaxis in wavenumber mode by using 'full' mode\n")
cmd = "aster = i2i(asu.data[i], from='lab1', to='aster', cm=1, full=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("plot lab and resampled ASTER spectra\n")
cmd = "pplot({asu.data[i], aster.data}, {'Lab Spectrum', 'ASTER'}, {1,5}, xaxis={asu.xaxis, aster.xaxis}, xlabel = 'Wavenumber', ylabel = 'Emissivity', key = '1400, .95')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("\n\n*** Example of how to do a bunch of instruments and plot vs wavelength\n\n")
printf("create aster x axis in wavelength mode (default is wavelength for broadband insts.)\n")
cmd = "xaster = make_band(inst='aster')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("create tims x axis in wavelength mode (default)\n")
cmd = "xtims = make_band(inst='tims')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("create themis x axis in wavelength mode (default)\n")
cmd = "xthemis = make_band(inst='themis')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("create lab1 x axis in wavelength mode\n")
cmd = "xlab1wave = make_band(inst='lab1', wave=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("create tes x axis in wavelength mode\n")
cmd = "xteswave = make_band(inst='tes', wave=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Convert lab spectrum to lab resolution with wavelength xaxis\n")
cmd = "lab1w = i2i(asu.data[i], from = 'lab1', to = 'lab1', wave=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to TES resolution with wavelength xaxis (default)\n")
cmd = "tes = i2i(asu.data[i], from = 'lab1', to = 'tes', wave=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to TIMS resolution with wavelength xaxis (default)\n")
cmd = "tims = i2i(asu.data[i], from = 'lab1', to = 'tims')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to TIMS resolution with wavelength xaxis (default)\n")
cmd = "themis = i2i(asu.data[i], from = 'lab1', to = 'themis')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to ASTER resolution with wavelength xaxis (default)\n")
cmd = "aster = i2i(asu.data[i], from = 'lab1', to = 'aster')"
echo(cmd)
eval(cmd)
pause("pausing")
printf("plot lab and resampled TES, THEMIS, ASTER, and TIMS spectra\n")
cmd = "pplot({lab1w, tes, themis, aster, tims}, {'Lab Spectrum', 'TES', 'THEMIS', 'ASTER', 'TIMS'}, {1,2,4,5,6}, xaxis={xlab1wave, xteswave, xthemis, xaster, xtims}, xlabel = 'Wavelength', ylabel = 'Emissivity', key = '12.5, .94', x1=6., x2=16.)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("\n\n**** Again - a simplier approach using 'full' mode in i2i ****\n\n")
printf("Resample lab spectrum to lab resolution with wavelength xaxis (default)\n")
cmd = "lab = i2i(asu.data[i], from = 'lab1', to = 'lab1', wave=1, full=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to TES resolution with wavelength xaxis (default)\n")
cmd = "tes = i2i(asu.data[i], from = 'lab1', to = 'tes', wave=1, full=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to THEMIS resolution with wavelength xaxis (default)\n")
cmd = "themis = i2i(asu.data[i], from = 'lab1', to = 'themis', wave=1, full=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to ASTER resolution with wavelength xaxis (default)\n")
cmd = "aster = i2i(asu.data[i], from = 'lab1', to = 'aster', wave=1, full=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("Resample lab spectrum to TIMS resolution with wavelength xaxis (default)\n")
cmd = "tims = i2i(asu.data[i], from = 'lab1', to = 'tims', wave=1, full=1)"
echo(cmd)
eval(cmd)
pause("pausing")
printf("plot lab and resampled TES, THEMIS, ASTER, and TIMS spectra\n")
cmd = "pplot({lab.data, tes.data, themis.data, aster.data, tims.data}, {'Lab Spectrum', 'TES', 'THEMIS', 'ASTER', 'TIMS'}, {1,2,4,5,6}, xaxis={lab.xaxis, tes.xaxis, themis.xaxis, aster.xaxis, tims.xaxis}, xlabel = 'Wavelength', ylabel = 'Emissivity', key = '12.5, .94', x1=6., x2=16.)"
echo(cmd)
eval(cmd)
pause("pausing")
}
define reverse_spectrum() {
if ($ARGC == 0) {
printf("Reverse the z dimension of an array\n")
printf("Usage: b = reverse_spectrum(a)\n")
printf(" p. christensen 12/08\n")
return(0)
}
spec = $1
nband = dim(spec)[3]
tmpspec = spec
for(i=1; i<=nband; i+=1) {
j = int((nband - (i-1)))
spec[,,j] = tmpspec[,,i]
}
return(spec)
}
define cm2micron() {
# convert wavenumber (input in cm) to micron
if ($ARGC == 0) {
printf("Convert array in cm-1 to microns\n")
printf("Usage: b = cm2micron(a)\n")
return(0)
}
# p. christensen 11/26/08
wave = 10000./$1
return(wave)
}
define r2t_hi(v){
if ($ARGC == 0) {
printf ("converts PRT resistance meassurments to temperature for bbhot \n")
printf ("(channels 3 and 4) using equations from ROSEMOUNT manual; this \n")
printf ("version is 1.5 that applies to PRTs with s/n 32147 and 31988 \n")
printf ("usage: r2t_hi(ch_3, ch_4) \n")
printf ("Written by Alice Baldridge - Thu Jun 15 13:41:16 MST 2006\n")
return(0)
}
if(HasValue(v)==0) {
v=1
}
#constants
A3=3.946065e-3
B3=-5.8755e-7
A4=3.947985e-3
B4=-5.8755e-7
Ro3=1001.01647
Ro4=1000.91666
#calculated values
W3=$1/Ro3
tp3 = ((A3^2-4*B3*(1-W3))^.5 - A3)/(2*B3)
r2t_hi_v1_temp3 = tp3+0.045*(tp3/100)*(tp3/100-1)*(tp3/419.58-1)*(tp3/630.74-1)
W4 = $2/Ro4
tp4 = ((A4^2-4*B4*(1-W4))^.5 - A4)/(2*B4)
r2t_hi_v1_temp4 = tp4+0.045*(tp4/100)*(tp4/100-1)*(tp4/419.58-1)*(tp4/630.74-1)
r2t_hi_v1_ave = (r2t_hi_v1_temp3+r2t_hi_v1_temp4)/2
if(v==1) printf ("Average temperature is \n")
return (r2t_hi_v1_ave)
}
define r2t_lo(v){
if ($ARGC == 0) {
printf ("converts PRT resistance meassurments to temperature for bbwarm \n")
printf ("(channels 1 and 2) using equations from ROSEMOUNT manual; this \n")
printf ("version is 1.5 that applies to PRTs with s/n 31985 and 31986 \n")
printf ("usage: r2t_hi(ch_1, ch_2) \n")
printf ("Written by Alice Baldridge - Thu Jun 15 13:41:16 MST 2006\n")
return(0)
}
if(HasValue(v)==0) {
v=1
}
# Constants
A1 = 3.946495e-3
B1 = -5.87549e-7
A2 = 3.947915e-3
B2 = -5.8755e-7
Ro1 = 1000.71666
Ro2 = 1001.59058
# Calculated values
W1 = $1/Ro1
tp1 = ((A1^2-4*B1*(1-W1))^.5 - A1)/(2*B1)
r2t_lo_temp1 = tp1+0.045*(tp1/100)*(tp1/100-1)*(tp1/419.58-1)*(tp1/630.74-1)
W2 = $2/Ro2
tp2 = ((A2^2-4*B2*(1-W2))^.5 - A2)/(2*B2)
r2t_lo_temp2 = tp2+0.045*(tp2/100)*(tp2/100-1)*(tp2/419.58-1)*(tp2/630.74-1)
r2t_lo_ave = (r2t_lo_temp1+r2t_lo_temp2)/2
if(v==1) printf("Average temperature is \n")
return(r2t_lo_ave)
}
define r2t_swri(v){
if ($ARGC < 4) {
printf("\nConverts PRT resistance measurements (ohms) to temperature (C) for BB #4 (SwRI)\n")
printf("using equations from the ROSEMOUNT manual \n\n")
printf("example: rt2_swri(ch1_lo,ch2_lo,ch1_hi,ch2_hi) \n\n")
printf("Osterloo/Hamilton 3/29/2010 \n\n")
return(NULL)
}
if(HasValue(v)==0) v=1
#Constants for BB #4 - PRT S/Ns 31991 and 31992 */
A1 = 3.947075e-3
B1 = -5.8755e-7
A2 = 3.946875e-3
B2 = -5.8755e-7
Ro1 = 1001.65018
Ro2 = 1000.40883
#Calculated values
W1 = $1/Ro1
tp1 = ((A1^2.0-4.0*B1*(1.0-W1))^0.5 - A1)/(2.0*B1)
r2t_lo_temp1 = tp1+0.045*(tp1/100.0)*(tp1/100.0-1.0)*(tp1/419.58-1.0)*(tp1/630.74-1.0)
W2 = $2/Ro2
tp2 = ((A2^2.0-4.0*B2*(1.0-W2))^.5 - A2)/(2.0*B2)
r2t_lo_temp2 = tp2+0.045*(tp2/100.0)*(tp2/100.0-1.0)*(tp2/419.58-1.0)*(tp2/630.74-1.0)
r2t_lo_avg = (r2t_lo_temp1+r2t_lo_temp2)/2.0
W3 = $3/Ro1
tp1_hi = ((A1^2.0-4.0*B1*(1.0-W3))^.5 - A1)/(2.0*B1)
r2t_hi_temp1 = tp1_hi+0.045*(tp1_hi/100.0)*(tp1_hi/100.0-1.0)*(tp1_hi/419.58-1.0)*(tp1_hi/630.74-1.0)
W4 = $4/Ro2
tp2_hi = ((A2^2.0-4.0*B2*(1.0-W4))^.5 - A2)/(2.0*B2)
r2t_hi_temp2 = tp2_hi+0.045*(tp2_hi/100.0)*(tp2_hi/100.0-1.0)*(tp2_hi/419.58-1.0)*(tp2_hi/630.74-1.0)
r2t_hi_avg = (r2t_hi_temp1+r2t_hi_temp2)/2.0
r2t_total = (r2t_lo_avg +r2t_hi_avg)/2.0
if(v!=0) printf ("Average temperature is %f\n",r2t_total)
return (r2t_total)
}
define emcal2(hotbb,warmbb,individual,ratw,dehyd,resample,noise_free,lab,max_emiss) {
if($ARGC==0) {
printf("\n")
printf ("emcal2() - Emissivity Calculation of Laboratory Spectra\n\n")
printf ("\tComputes an instument response function from measured emission, \n")
printf ("\tspectra of hot and warm blackbody and known blackbody Temperature. \n")
printf ("\tThen calculates emissivity using previously determined\n")
printf ("\tresponse function, hot blackbody spectrum, sample spectrum,\n")
printf ("\tand known environment temperature.\n\n")
printf ("Usage:\n")
printf ("\t$1=directory containing spectra\n")
printf ("\t$2=output file name\n")
printf ("\t$3=warm BB resistance values (Ch 101//Ch102 as 1313.3//1315.6(asu) or Ch1lo//Ch2lo//Ch1hi//Ch2hi(swri))\n")
printf ("\t$4=hot BB resistance values (Ch 103//Ch104 as 1314.3//1314.6(asu) or Ch1lo//Ch2lo//Ch1hi//Ch2hi(swri))\n")
printf ("Or\n")
printf ("\t$3=warm BB temperture (from r2t_lo() as 70.403)\n")
printf ("\t$4=hot BB temperture (from r2t_hi() as 100.334)\n\n")
printf ("Alternatively, if using the ASU LabView data acquisition software\n")
printf ("\t$1=the full path to the data qube written by the labview interface \n")
printf ("\t$2=the output file name (which will be written into the data qube directory\n")
printf ("Optional Arguments:\n")
printf ("\t$5=\tlist of environmental chamber temperatures\n")
printf ("\t\tmust have the same # of entries as spectra in $1\n")
printf ("\tlab=\tlocation where data was acquired (PRT calibration values are lab specific) (\"swri\" or Default=\"asu\")\n")
printf ("\tnoise_free=compute a \"noise free\" instrument response function\n\t\tand instrument radiance (Default=1)\n")
printf ("\thotbb=\thot blackbody partial filename (Default search = \"bbh\")\n")
printf ("\twarmbb=\twarm blackbody partial filename (Default search = \"bbw\")\n")
printf ("\tdehyd=\tadd dehydrated spectra to the output\n\t\tstructure (Uses dehyd2(), Default=0)\n")
printf ("\tindivdual=write the spectra out and individual text files with\n\t\tno labels for spectral library (Default=0)\n")
printf ("\tresample=resample all spectra to hot blackbody xaxis (Default=0) \n")
printf ("\tmax_emiss=optional maximum emissivity (e.g. 0.98) (default = 1.)\n")
printf ("\tratw=\tset up default options for RATW calibration (Default=0, deyhd=1,\n\t\tindiv=1 (plus modified label), resample=1) \n\n")
printf ("NOTE: This function uses irf_ri(), cal_rad(), and emissivity() to do calibration\n\n")
printf ("Written by C.Edwards. 3/18/2011\n\n")
return(null)
}
raw={}
out={}
dir=$1
outname=$2
if(HasValue(individual)==0) individual=0
if(HasValue(dehyd)==0) dehyd=0
if(HasValue(ratw)==0) ratw=0
if(HasValue(resample)==0) resample=0
if(HasValue(noise_free)==0) noise_free=1
if(HasValue(lab)==0) lab="asu"
if(HasValue(max_emiss)==0) max_emiss=1.
if(ratw==1) {
dehyd=0
lab="asu"
noise_free=1
individual=1
resample=1
}
#check to make sure the output file is a string and add a .hdf if it doesn't exist
if(type(outname)!="STRING"){
printf("Please use a string for output filename\n")
return(NULL)
}
if(strstr(outname, ".hdf")== 0) {
outname = outname+".hdf"
}
raw.hdfname=outname
raw.inputdir=dir
if(type($1)=="STRING" && $ARGC==2) {
verbose=0
input=load_pds(dir)
verbose=3
if(HasValue(get_struct(input,"spectral_data_table"))==0) {
printf("\nThe input file you specify is not a valid spectral qube\n")
printf("Please check your path and verify the input file name\n\n")
return(null)
}
raw.hotbbfile=dir
raw.warmbbfile=dir
raw.inputdir=dirname(dir)+"/"
spec_data_size=int(input.spectral_data_table.data.spec_data_size[,1]-2)
raw.data=create(1,1,spec_data_size,format=float,org=bsq,start=0,step=0)
raw.hotbb=create(1,1,spec_data_size,format=float,org=bsq,start=0,step=0)
raw.warmbb=create(1,1,spec_data_size,format=float,org=bsq,start=0,step=0)
envtemp=create(1,1,1,format=float,org=bsq,start=0,step=0)
raw.insttemp=create(1,1,1,format=float,org=bsq,start=0,step=0)
raw.speclist=text(0)
for(i=1;i<=length(input.spectral_data_table.data.acq_type);i+=1) {
if(input.spectral_data_table.data.acq_type[,i]=="S") {
raw.data=cat(raw.data,float(translate(input.spectral_data_table.data.spec_data[2:spec_data_size+1,i],x,z)),axis=x)
envtemp=cat(envtemp,float(input.spectral_data_table.data.env_temp[,i]),axis=x)
raw.insttemp=cat(raw.insttemp,float(input.spectral_data_table.data.inst_temp[,i]),axis=x)
raw.speclist=cat(raw.speclist,strsub(input.spectral_data_table.data.label[,i]," ",""),axis=y)
}
if(input.spectral_data_table.data.acq_type[,i]=="W") {
raw.warmbb=float(translate(input.spectral_data_table.data.spec_data[2:spec_data_size+1,i],x,z))
warmres=input.spectral_data_table.data.bb_res_1[,i]//input.spectral_data_table.data.bb_res_2[,i]
}
if(input.spectral_data_table.data.acq_type[,i]=="H") {
raw.hotbb=float(translate(input.spectral_data_table.data.spec_data[2:spec_data_size+1,i],x,z))
hotres=input.spectral_data_table.data.bb_res_1[,i]//input.spectral_data_table.data.bb_res_2[,i]
}
}
out.xaxis=translate(input.spectral_xaxis[2:spec_data_size+1],x,z)
raw.data=raw.data[2:]
envtemp=envtemp[2:]
raw.insttemp=raw.insttemp[2:]
raw.input=input
lab="asu"
} else {
warmres=$3
hotres=$4
#if we have a 5th argument it is the enviornment temperature array in x...
if($ARGC>4) {
if(dim($5)[2]>1 && dim($5)[1]==1) {
envtemp=float(translate($5,x,y))
} else if(dim($5)[2]==1 && dim($5)[1]>1){
envtemp=$5
} else {
printf("\nInvalid environmental temperature array supplied\n")
printf("Please re-enter the environmental array or remove the option\n")
printf("to enter the temperatures manually.\n\n")
return(null)
}
}
#set up the default variables
if(HasValue(hotbb)==0) {
hotbb="bbh"
} else {
hotbb=basename(hotbb,".CSV")
}
if(HasValue(warmbb)==0) {
warmbb="bbw"
} else {
warmbb=basename(warmbb,".CSV")
}
#define the lab where spectra were acquired
out.lab=lab
#check to make sure the input directory is a string
#also check that it exists
#test to make sure the directory ends with a /
if(type(dir)!="STRING"){
printf("Please use a string for the input directory\n")
return(NULL)
}
if(fexists(dir)==0) {
printf("Directory does not exist: %s\n",dir)
printf("Please enter a valid directory\n")
return(NULL)
}
if(basename(dir)==".") dir=$PWD
if(dir[length(dir)]!="/") dir=dir+"/"
#find the apporpriate files
list=syscall("ls "+dir+"*.CSV")
raw.speclist=text()
for(i=1;i<=length(list);i+=1) {
if(strstr(list[,i],warmbb)==0 && strstr(list[,i],hotbb)==0) {
raw.speclist=cat(raw.speclist,list[,i],axis=y)
} else if (strstr(list[,i],warmbb)!=0 ) {
raw.warmbbfile=list[,i]
} else if (strstr(list[,i],hotbb)!=0 ) {
raw.hotbbfile=list[,i]
}
}
#test to see that we found the blackbodies
if(HasValue(get_struct(raw,"hotbbfile"))==0 || HasValue(get_struct(raw,"warmbbfile"))==0) {
printf("\nCould not file the hot or warm blackbody file\n")
printf("Check then names and try again\n\n")
return(null)
}
verbose=0
#read in the black bodies and store an xaxis from the hot blackbody...
hotbb=translate(ascii(raw.hotbbfile,format=float,delim=","),y,z)
raw.hotbb=hotbb[2]
raw.warmbb=translate(ascii(raw.warmbbfile,format=float,delim=",")[2],y,z)
out.xaxis=hotbb[1]
#load the data, resampling if needed and user specified
raw.data=create(length(raw.speclist),1,dim(out.xaxis)[3],format=float,org=bsq)
for(i=1;i<=length(raw.speclist);i+=1) {
tmpdata=translate(ascii(raw.speclist[,i],delim=",",format=float),y,z)
if(resample==0 && dim(tmpdata)[3]!=dim(out.xaxis)[3]) {
printf("\nYour spectral data and xaxis do not match for at least one sample\n")
printf("Please rerun with resample=1\n\n")
verbose=3
return(null)
} else if(resample==1 && dim(tmpdata)[3]!=dim(out.xaxis)[3]) {
raw.data[i]=resample(oldy=tmpdata[2],oldx=tmpdata[1],newx=out.xaxis).data
} else {
raw.data[i]=tmpdata[2]
}
}
verbose=3
}
#check to see if we have resistances or temperatures
if(sum(dim(warmres))==3 && sum(dim(hotres))==3) {
raw.warmbbt=warmres
raw.hotbbt=hotres
} else {
if(lab=="asu") {
#check the dimensions of the Warm BB resistance values
if(isnum(warmres)==0 || dim(warmres)[1]!=2) {
printf("Please enter resistance values for the WarmBB as ch1//ch2 as $3\n")
return(NULL)
}
#check the dimensions of the Hot BB resistance values
if(isnum(hotres)==0 || dim(hotres)[1]!=2) {
printf("Please enter resistance values for the HotBB as ch3//ch4 as $4\n")
return(NULL)
}
#add these to the output struct
raw.ch1res=float(warmres[1])
raw.ch2res=float(warmres[2])
raw.ch3res=float(hotres[1])
raw.ch4res=float(hotres[2])
#convert the resistance values to temperatures
raw.warmbbt=r2t_lo(raw.ch1res,raw.ch2res,v=0)
raw.hotbbt=r2t_hi(raw.ch3res,raw.ch4res,v=0)
} else if (lab=="swri") {
#check the dimensions of the Warm BB resistance values
if(isnum(warmres)==0 || dim(warmres)[1]!=4) {
printf("Please enter resistance values for the WarmBB as ch1_lo//ch2_lo//ch1_hi//ch2_hi as $3\n")
return(NULL)
}
#check the dimensions of the Hot BB resistance values
if(isnum(hotres)==0 || dim(hotres)[1]!=4) {
printf("Please enter resistance values for the HotBB as ch1_lo//ch2_lo//ch1_hi//ch2_hi as $4\n")
return(NULL)
}
#add these to the output struct
raw.warm_ch1res_lo=float(warmres[1])
raw.warm_ch2res_lo=float(warmres[2])
raw.warm_ch1res_hi=float(warmres[3])
raw.warm_ch2res_hi=float(warmres[4])
raw.hot_ch1res_lo=float(hotres[1])
raw.hot_ch2res_lo=float(hotres[2])
raw.hot_ch1res_hi=float(hotres[3])
raw.hot_ch2res_hi=float(hotres[4])
#convert the resistance values to temperatures
raw.warmbbt=r2t_swri(raw.warm_ch1res_lo,raw.warm_ch2res_lo,raw.warm_ch1res_hi,raw.warm_ch2res_hi,v=0)
raw.hotbbt=r2t_swri(raw.hot_ch1res_lo,raw.hot_ch2res_lo,raw.hot_ch1res_hi,raw.hot_ch2res_hi,v=0)
} else {
printf("No valid lab location was specified\n")
printf("Pleas re-enter either \"asu\" or \"swri\"\n\n")
return(null)
}
}
#now we see if the user needs to input data or if it was supplied beforehand
userinput=0
if(HasValue(envtemp)) {
if(dim(envtemp)[1]!=dim(raw.data)[1]) {
printf("\nThe number of enviornement temperatures does not equal the number of spectra\n")
printf("Either remove the environment argument to enter manually, or fix the input\n")
printf("# of Spectra:\t\t%i\n# of Env Temperatures: \t%i\n\n",length(raw.speclist),dim(envtemp)[1])
return(null)
} else {
raw.envtemp=envtemp
}
} else {
userinput=1
}
#not supplied beforehand, get ready for some manual entry!!
if(userinput==1) {
raw.envtemp=create(length(raw.speclist),1,1,format=float,org=bsq,start=0,step=0)
printf("\nBeginning manual enviornmental temperature entry\n")
printf("Enter \"q\" to quit at any time\n")
for(i=1;i<=length(raw.speclist);i+=1) {
val=pause(sprintf("\nEnter the Chamber Temperature (Ch105 in degrees C) for\nSample: %s\n",basename(raw.speclist[,i],".CSV")))
if(val=="\n") {
i-=1
} else if (val[1,1]=="q" || val[1,1]=="Q") {
printf("\nExiting...\n\n")
return(null)
} else {
raw.envtemp[i]=atof(val)
}
}
}
out.raw=raw
if(noise_free!=0) {
#use irf_ri() to compute the noisefree instrument response function and instrument contributed radiance
out.irf_ri=irf_ri(bb1=raw.hotbb,t1=ctok(raw.hotbbt),bb2=raw.warmbb,t2=ctok(raw.warmbbt),xaxis=out.xaxis,noise_free=1)
} else {
#use irf_ri() to compute the standard instrument response function and instrument contributed radiance
out.irf_ri=irf_ri(bb1=raw.hotbb,t1=ctok(raw.hotbbt),bb2=raw.warmbb,t2=ctok(raw.warmbbt),xaxis=out.xaxis)
}
#use cal_rad() to compute calibrated radaince from a spectrum
out.cal_rad=cal_rad(spec=raw.data,ri=out.irf_ri.ri,irf=out.irf_ri.irf,xaxis=out.xaxis)
#use emissivity() to compute emissivity from calibrated radiance
out.emissivity=emissivity(out.cal_rad.radiance, downwelling_t=ctok(raw.envtemp),xaxis=out.xaxis,max_emiss=max_emiss)
#set up some convineient fields for the output structure
out.radiance=out.cal_rad.radiance
out.sample_temperature=ktoc(out.emissivity.target_temperature)
out.sample_t_channel=out.emissivity.pos_t1
out.sample_t_channel[where out.emissivity.target_temperature<=out.emissivity.threshold_t_warm]=out.emissivity.pos_t2
out.sample_t_wavenumber=translate(subset(out.xaxis,z=out.sample_t_channel),z,x)
out.data=out.emissivity.emissivity
#apply dehydration algorithm if the user specifies it
if(dehyd==1) {
dehyd={}
dehyd.original_data=out.data
for(i=1;i<=dim(out.data)[1];i+=1) {
out.data[i]=dehyd2(dehyd.original_data[i],out.xaxis)
}
out.dehyd=dehyd
}
#if ratw is selected then add the emcal emissivity of....to the label
if(ratw==1) {
out.label="emcal emissivity of "+basename(raw.speclist,".CSV")
} else {
out.label=basename(raw.speclist,".CSV")
}
#if individual is selected then write out vm compatible text files
if(individual==1) {
for(i=1;i<=length(out.label);i+=1) {
indiv={}
indiv.data=float(out.data[i])
indiv.xaxis=float(out.xaxis)
indiv.label=out.label[,i]
write_vm(indiv,raw.inputdir+basename(raw.speclist[,i],".CSV")+".txt")
}
}
if(ratw==1) {
out.label=basename(raw.speclist,".CSV")
}
#output the data
write(out,raw.inputdir+raw.hdfname,hdf,force=1)
return(out)
}
define emcal(cold, indiv, libspec, field_spec, resample, ratw, bbhot, bbwarm, dehyd) {
#added support for other names than standard bb names - c.edwards 1-11-08
#added checking for .hdf output suffix - c.edwards 1-11-08
#added dehyd2 support - c.edwards 1-11-08
#fixed spacing/verbosity for ease of viewing - c.edwards 1-11-08
#fixed indiv to write out individual text files- c.edwards 1-11-08
if($ARGC <4) {
printf ("emcal() - Emissivity Calculation of Laboratory Spectra\n\n")
printf (" Computes an instument response function from measured emission, \n")
printf (" spectra of hot and warm blackbody and known blackbody T. \n")
printf (" Then calculates emissivity using previously calculated\n")
printf (" response function, hot bb spectrum, hot sample spectrum,\n")
printf (" and known environment temperature.\n\n")
printf (" NOTE:REQUIRES thm module \n\n")
printf ("usage: $1=directory containing spectra \n")
printf (" $2=hot bb temp \n")
printf (" $3=warm bb temp \n")
printf (" $4=output file name\n")
printf ("\t Optional Arguments\n")
printf (" bbhot=hot black body filename (Default=\"bbhot.CSV\")\n")
printf (" bbwarm=warm black body filename (Default=\"bbwarm.CSV\")\n")
printf (" dehyd=add dehydrated spectra to the output structure (Uses dehyd2()) Default=\"no\"\n")
printf (" indiv=write the spectra out and individual text files with no labels for spectral library Default=\"no\"\n")
printf (" libspec=prompt for input data for spectral library\n")
printf (" resample=resample spectra to xaxis incase OMEGA saved and extra few \n")
printf (" cold=some or all files measured below chamber temperature\n")
printf (" field_spec=calibrate D&P spectra \n")
printf (" ratw=RATW files \n")
printf (" temp=prompt for temperature when sample temperature is known in degrees C\n")
printf (" wn=prompt for wavenumber for sample determination\n\n")
printf (" Written by Alice Baldridge - Thu Jun 15 13:41:16 MST 2006\n\n")
return(NULL)
}
directory = $1
directory=directory+"/"
if(fexists(directory)==0) {
printf("Directory does not exist: %s\n",directory)
printf("Please enter a valid directory\n")
return(NULL)
}
hot_bb_t=$2
if(isnum(hot_bb_t)==0) {
printf("Please enter a number for the Hot Black Body Temperature\n")
return(NULL)
}
warm_bb_t=$3
if(isnum(warm_bb_t)==0) {
printf("Please enter a number for the Warm Black Body Temperature\n")
return(NULL)
}
out=$4
if(type(out)!="STRING"){
printf("please use a string for output filename\n")
return(NULL)
}
pd = strstr(out, ".hdf")
if(pd == 0) {
out = out+".hdf"
}
deh="no"
if (HasValue(dehyd)) deh=dehyd
rat="no"
if (HasValue(ratw)) rat=ratw
ind="no"
if (HasValue(indiv)) ind=indiv
brr="no"
if (HasValue(cold)) brr=cold
lib="no"
if (HasValue(libspec)) lib=libspec
field="no"
if (HasValue(field_spec)) field=field_spec
resamp="no"
if (HasValue(resample)) resamp=resample
if(field=="yes"){
batch_extract(directory)
list =basename( syscall("ls -t "+directory+"*.txt"))
bbhot_name = directory + "WARMBB.WBB.txt"
bbwarm_name = directory + "COLDBB.CBB.txt"
gold_name = directory + "GOLD.DWR.txt"
}else{
list =basename( syscall("ls -t "+directory+"*.CSV"))
bbhot_name = directory + "bbhot.CSV"
bbwarm_name = directory + "bbwarm.CSV"
}
#added support for other names than standard bb names - c.edwards 1-11-08
if(HasValue(bbhot)) bbhot_name=directory+basename(bbhot)
if(HasValue(bbwarm)) bbwarm_name=directory+basename(bbwarm)
if(fexists(bbhot_name)==0 || fexists(bbwarm_name)==0) {
printf("Hot or Warm black body name is incorrect\n")
return(null)
}
len = length(list)
####################################
# read in the black bodies and xaxis
####################################
verbose=0
# structure containing the spectra
spec_struct = {}
# read in bbhot
if(field=="yes"){
tmp=ascii(bbhot_name, format=float)
}else{
tmp = ascii(bbhot_name,format=float,delim=",")
}
# add Xaxis to the spec_struct
spec_struct.xaxis = translate(tmp[1],from=y,to=z)
# add bbhot to spec_struct
spec_struct.hbb = translate(tmp[2],from=y,to=z)
spec_struct.hotbb_temperature=hot_bb_t
# read in bbwarm
if(field=="yes"){
tmp=ascii(bbwarm_name, format=float)
}else{
tmp = ascii(bbwarm_name,format=float,delim=",")
}
# add bbwarm to spect_struct
spec_struct.wbb = translate(tmp[2],from=y,to=z)
spec_struct.warmbb_temperature=warm_bb_t
# read in the gold plate
if(field=="yes"){
tmp=ascii(gold_name, format=float)
spec_struct.gold= translate(tmp[2],from=y,to=z)
}
##############################
# create the response function
##############################
spec_struct.response=(spec_struct.hbb-spec_struct.wbb)/(bbr(spec_struct.xaxis,ctok(hot_bb_t))-bbr(spec_struct.xaxis,ctok(warm_bb_t)))
spec_struct.response=float(spec_struct.response)
xplot(spec_struct.response[1,,],Xaxis=spec_struct.xaxis)
plim(0,4000,0,0)
nullim()
lines()
pause("\nDisplaying Response Function; Press \"Enter\" to continue")
#################################
# fill in the raw data and labels
#################################
spec_names = ""
spectra = spec_struct.xaxis*0
for(i=1;i<=len;i+=1) {
filename = list[,i]
if(field=="yes"){
if(filename !=basename(bbwarm_name) && filename != basename(bbhot_name) && filename != "GOLD.DWR.txt" && filename != "*.asc" && filename != "bbhot.CSV" && filename != "bbwarm.CSV"){
data = translate(ascii(directory+filename,format=float),from=y,to=z)
data = data[2]
extension = strstr(filename,".txt")
filename = filename[:extension-1]
# add to the data block
spec_names = cat(spec_names,filename,axis=y)
spectra = cat(spectra,data,axis=y)
}
}else{
if(filename !=basename(bbhot_name) && filename != basename(bbwarm_name) && filename != "bbhot.CSV" && filename != "bbwarm.CSV") {
data = translate(ascii(directory+filename,format=float,delim=","),from=y,to=z)
xtemp=data[1]
data = data[2]
if (resamp=="yes"){
data = resample(data, xtemp, spec_struct.xaxis).data
}
extension = strstr(filename,".CSV")
filename = filename[:extension-1]
# add to the data block
spec_names = cat(spec_names,filename,axis=y)
spectra = cat(spectra,data,axis=y)
}
}
}
verbose=3
# trim off the blank space of our spectra and spec_names
spec_names = spec_names[,2:]
spectra = spectra[,2:]
# paste these motherfuckers into the structure
spec_struct.raw = spectra
spec_struct.sample_name = spec_names
###################################
# radiance and emissivity calculation
###################################
emissivity=spec_struct.xaxis*0
radiance=spec_struct.xaxis*0
btemp=spec_struct.xaxis*0
dm = dim(spec_struct.raw)[2]
sample_temperature=""
chamber_temperature=""
sample_id=""
quality=""
category=""
latitude=""
longitude=""
location=""
locality=""
particle_size=""
citation=""
type=""
mineral_number=""
chemical_formula=""
subgroup=""
#calculate the gold plate radiance
if(field=="yes"){
gold_T=pause("Enter the gold plate temperature\n")
gold_T=atof(gold_T[,1])
gold_energy=bbr(spec_struct.xaxis, ctok(gold_T))
B_inst=spec_struct.hbb-(spec_struct.response*(bbr(spec_struct.xaxis,ctok(hot_bb_t))))
goldrad = (spec_struct.gold - B_inst)/spec_struct.response
factor=0.04
Ldwr= (goldrad-factor*gold_energy)/(1.0-factor)
}
for(k=1;k<=dm;k+=1) {
name=spec_struct.sample_name[,k]
#determine sample radiance
if(field=="yes"){
B_inst=spec_struct.hbb-(spec_struct.response*(bbr(spec_struct.xaxis,ctok(hot_bb_t))))
samp_T=pause("Enter the sample temperature\n")
samp_T=atof(samp_T[,1])
samp_energy=bbr(spec_struct.xaxis, ctok(samp_T))
rad = (spec_struct.raw[,k] - B_inst)/spec_struct.response
rad = float(rad)
}else{
B_inst=(bbr(spec_struct.xaxis,ctok(hot_bb_t))-(spec_struct.hbb/spec_struct.response))
inst_T=avg(btemp(spec_struct.xaxis[,,208:520],B_inst[,,208:520]))
printf("\nInstrument T is: %.2f\n", ktoc(inst_T))
energy=(bbr(spec_struct.xaxis,inst_T))
if(brr=="yes"){
opt_c=pause("Was "+name+" measured cold? \n")
if(opt_c[1,1] == "y") {
spec_struct.raw[,k]= (-1)*spec_struct.raw[,k]
}
}
rad = ((spec_struct.raw[,k]/spec_struct.response) + energy)
rad = float(rad)
}
#########################
#emissivity calculations
#########################
#determine sample maximum brightnes temperature
B_temp=btemp(spec_struct.xaxis,rad)
B_temp=float(B_temp)
if(field=="yes"){
kern=clone(1,1,20,1)
samp_T=max(thm.convolve(B_temp[1,,],kern)[,,250:500])
n=maxchan(thm.convolve(B_temp[1,,],kern)[,,250:500])+250
samp_T_new=ktoc(samp_T)
emcal_T_wavenumber=spec_struct.xaxis[,,n]
sample_temp=sprintf("%.2f\n",samp_T_new)
samp_T_new=ctok(samp_T_new)
sample_temp=sample_temp[1:(strlen(sample_temp)-1)]
sample_temperature=cat(sample_temperature, sample_temp, axis=y)
}else{
samp_T_new=ctok(30)
kern=clone(1.0,1,1,26)
samp_T=max(thm.convolve(B_temp[1,,],kern)[,,50:831])
if(brr=="yes"){
if(opt_c[1,1] == "y"){
samp_T_new=min(thm.convolve(B_temp[1,,],kern)[,,105:831])
m=minchan(thm.convolve(B_temp[1,,],kern)[,,105:831])+105
emcal_T_wavenumber=spec_struct.xaxis[,,m]
}else{
samp_T_new=samp_T
n=maxchan(thm.convolve(B_temp[1,,],kern)[,,105:831])+105
emcal_T_wavenumber=spec_struct.xaxis[,,n]
}
}else{
samp_T_new=samp_T
n=maxchan(thm.convolve(B_temp[1,,],kern)[,,50:831])+50
emcal_T_wavenumber=spec_struct.xaxis[,,n]
}
printf("Sample T for \""+name+"\" is: %.2f\n",ktoc(samp_T_new))
sample_t=ktoc(samp_T_new)
sample_temp=sprintf("%.2f\n",sample_t)
sample_temp=sample_temp[1:(strlen(sample_temp)-1)]
sample_temperature=cat(sample_temperature, sample_temp, axis=y)
}
printf("Sample T for sample \""+name+"\" is determined at: %.2f cm^-1\n",emcal_T_wavenumber)
B_samp = (bbr(spec_struct.xaxis,samp_T_new))
#determine emissivity
if(field=="yes"){
emis=(rad - Ldwr)/(B_samp-Ldwr)
emis = float(emis)
chamber_temperature="na"
}else{
#determine environmental radiance
env_T=pause("What is the chamber T for \""+name+"\"?\n")
chamber_t=env_T[1:(strlen(env_T)-1)]
if(env_T[1,1]=="q"){
return(0)
}else{
env_T=atof(env_T[,1])
env_rad=(bbr(spec_struct.xaxis,ctok(env_T)))
}
chamber_temperature=cat(chamber_temperature, chamber_t, axis=y)
emis=(((spec_struct.raw[,k]/spec_struct.response)-env_rad+energy)/(B_samp - env_rad))
emis = float(emis)
#water removal for RATW
if (rat=="yes"){
emis=dehyd2(emis,spec_struct.xaxis)
}
}
#plot
xplot(emis[1,,],Xaxis=spec_struct.xaxis)
lines()
plim(2000,200, 0.0, 1.05)
radiance=cat(radiance, rad, axis=y)
btemp=cat(btemp, B_temp, axis=y)
emissivity=cat(emissivity,emis,axis=y)
#####################################################
#enter the spec lib options for the individual files
#####################################################
if(lib=="yes"){
sampid=pause("Enter the sample id for "+name+" \n")
if(sampid == "\n"){
sampleid=""
}else{
sampleid = sampid[1:(strlen(sampid)-1)]
}
sample_id=cat(sample_id,sampleid, axis=y)
gory=pause("Enter the sample deconvolution category for "+name+" \n")
if (gory == "\n"){
egory = ""
}else{
egory = gory[1:(strlen(gory)-1)]
}
category=cat(category, egory, axis=y)
qual=pause("Enter the quality flag for "+name+": 1=Well characterized, pure, 2=Well characterized, minor impurities, 3=Well characterized, major impurities, 4=Not well characterizied\n")
if(qual[1,1]=="\n"){
qu="[4] Not well characterized"
}else{
if(qual[1,1]=="1"){
qu="[1] Well characterized, pure"
}else{
if(qual[1,1]=="2"){
qu="[2] Well characterized, minor impurities"
}else{
if(qual[1,1]=="3"){
qu="[3] Well characterized, major impurities"
}else{
if(qual[1,1]=="4"){
qu="[4] Not well characterized"
}}}}}
quality=cat(quality,qu,axis=y)
lat=pause("Enter the sample latitude for "+name+" \n")
if(lat=="\n"){
la=""
}else{
la=lat[1:(strlen(lat)-1)]
}
latitude=cat(latitude, la, axis=y)
long=pause("Enter the sample longitude for "+name+" \n")
if(long=="\n"){
lo=""
}else{
lo=long[1:(strlen(long)-1)]
}
longitude=cat(longitude, lo, axis=y)
loc=pause("Enter the sample location for "+name+" \n")
if(loc=="\n"){
tion=""
}else{
tion=loc[1:(strlen(loc)-1)]
}
location=cat(location, tion, axis=y)
local=pause("Enter the sample collection locality for "+name+" \n")
if(local=="\n"){
lity=""
}else{
lity=local[1:(strlen(local)-1)]
}
locality=cat(locality, lity, axis=y)
part=pause("Enter the particle size for "+name+" 1=Fine particulate, 2=Coarse particulate, 3=Pebble,4=Handsample\n")
if(part[1,1]=="\n"){
particlesize=""
}else{
if(part[1,1]=="1"){
particlesize="Fine particulate"
}else{
if(part[1,1]=="2"){
particlesize="Coarse particulate"
}else{
if(part[1,1]=="3"){
particlesize="Pebble"
}else{
if(part[1,1]=="4"){
particlesize="Handsample"
}}}}}
particle_size=cat(particle_size, particlesize, axis=y)
cit=pause("Enter the sample citation for "+name+" \n")
if(cit=="\n"){
cite=""
}else{
cite=cit[1:(strlen(cit)-1)]
}
citation=cat(citation, cite, axis=y)
typ=pause("Enter the sample type for "+name+" 1=Mineral/Phase, 2=Rock, 3=Soil/Regolith, 4=Meteorite/Lunar, 5=Other\n")
if(typ[1,1]=="\n"){
ty=""
}else{
if(typ[1,1]=="1"){
ty="Mineral/Phase"
printf("Enter ther dana mineral group for "+name+" \n")
printf("1=Native Elements, 2=Sulfides, 3=Oxides and Hydroxides, \n")
printf("4=Halides, 5=Carbonates, Borates and Nirtates, 6=Sulfates, \n")
printf("Chromates and Molydates, 7=Phosphates, Arsenates and Vanadates,\n")
printf("8=Organic Materials 9=Nesosilicates, 10=Sorosilicates, \n")
printf("11=Cyclosilicates, 12=Inosilicates, 13=Phyllosilicates, \n")
sub=pause("14=Tectosilicates, 15=Other Silicates, 16=Other \n")
if(sub[1,1]=="\n"){
sg=""
}else{
if(sub[1,1]=="1"){
sg="I - Native Elements"
}
if(sub[1:(strlen(sub)-1)]=="2"){
sg="II - Sulfides"
}
if(sub[1:(strlen(sub)-1)]=="3"){
sg="III - Oxides and Hydroxides"
}
if(sub[1:(strlen(sub)-1)]=="4"){
sg="IV - Halides"
}
if(sub[1:(strlen(sub)-1)]=="5"){
sg="V - Carbonates, Nitrates, Borates"
}
if(sub[1:(strlen(sub)-1)]=="6"){
sg="VI - Sulfates, Chromates, Molybdates"
}
if(sub[1:(strlen(sub)-1)]=="7"){
sg="VII - Phosphates, Arsenates, Vanadates"
}
if(sub[1:(strlen(sub)-1)]=="8"){
sg="IX - Organic Minerals"
}
if(sub[1:(strlen(sub)-1)]=="9"){
sg="VIII - Silicates - Nesosilicates (e.g. Olivine, Garnet)"
}
if(sub[1:(strlen(sub)-1)]=="10"){
sg="VIII - Silicates - Sorosilicates (e.g. Epidote)"
}
if(sub[1:(strlen(sub)-1)]=="11"){
sg="VIII - Silicates - Cyclosilicates (e.g. Tourmaline, Beryl)"
}
if(sub[1:(strlen(sub)-1)]=="12"){
sg="VIII - Silicates - Inosilicates (e.g. Pyroxene, Amphibole)"
}
if(sub[1:(strlen(sub)-1)]=="13"){
sg="VIII - Silicates - Phyllosilicates (e.g. Micas, Clays)"
}
if(sub[1:(strlen(sub)-1)]=="14"){
sg="VIII - Silicates - Tektosilicates (e.g. Quartz, Feldspars, Zeolites)"
}
if(sub[1:(strlen(sub)-1)]=="15"){
sg="VIII - Silicates - Other Silicates"
}
if(sub[1:(strlen(sub)-1)]=="16"){
sg="Other"
}}
minnum=pause("Enter the dana mineral number for "+name+" \n")
if(minnum=="\n"){
num=""
}else{
num=minnum[1:(strlen(minnum)-1)]
}
form=pause("Enter the chemical formula for "+name+" \n")
if(form=="\n"){
fo=""
}else{
fo=form[1:(strlen(form)-1)]
}
}
if(typ[1,1]=="2"){
ty="Rock"
num=""
fo=""
printf("Enter the rock subgroup for "+name+" \n")
printf("1=Igneous, 2=Sedimentary, 3=Metamorphic, \n")
sub=pause("4=Other\n")
if(sub[1,1]=="\n"){
sg=""
}else{
if(sub[1,1]=="1"){
sg="Igneous"
}else{
if(sub[1,1]=="2"){
sg="Sedimentary"
}else{
if(sub[1,1]=="3"){
sg="Metamorphic"
}else{
if(sub[1,1]=="4"){
sg="Other"
}}}}}
}
if(typ[1,1]=="3"){
ty="Soil/Regolith"
num=""
fo=""
sub=pause("Enter the soil/regolith subgroup for "+name+" \n")
if(sub=="\n"){
sg=""
}else{
sg=sub[1:(strlen(sub)-1)]
}
}
if(typ[1,1]=="4"){
ty="Meteorite/Lunar"
num=""
fo=""
printf("Enter the meteorite/lunar subgroup for "+name+" \n")
printf("1=Chondrite, 2=Primitive achondrite, 3=Asteroid/Lunar/Martian achondrite, \n")
sub=pause("4=Iron/stoney-iron achondrite \n")
if(sub[1,1]=="\n"){
sg=""
}else{
if(sub[1,1]=="1"){
sg="Chondrite"
}
if(sub[1,1]=="2"){
sg="Primitive achondrite"
}
if(sub[1,1]=="3"){
sg="Asteroid/Lunar/Martian achondrite"
}
if(sub[1,1]=="4"){
sg="Iron/Stoney-iron achondrite"
}}
}
if(typ[1,1]=="5"){
ty="Other"
}
}
type=cat(type, ty, axis=y)
mineral_number=cat(mineral_number, num, axis=y)
subgroup=cat(subgroup, sg, axis=y)
chemical_formula=cat(chemical_formula, fo, axis=y)
#close the lib=yes loop
}
#close the for loop
}
###############################
#Put everything in the structure
###############################
sample_temperature=sample_temperature[,2:]
spec_struct.sample_temperature=sample_temperature
chamber_temperature=chamber_temperature[,2:]
spec_struct.chamber_temperature=chamber_temperature
radiance= radiance[,2:]
spec_struct.radiance=radiance
btemp=btemp[,2:]
spec_struct.btemp=btemp
emissivity = emissivity[,2:]
spec_struct.data = emissivity
if (lib=="yes"){
sample_id=sample_id[,2:]
spec_struct.sample_id=sample_id
category=category[,2:]
spec_struct.category=category
quality = quality[,2:]
spec_struct.quality=quality
latitude=latitude[,2:]
spec_struct.latitude=latitude
longitude=longitude[,2:]
spec_struct.longitude=longitude
location=location[,2:]
spec_struct.sample_location=location
locality=locality[,2:]
spec_struct.collection_locality=locality
particle_size=particle_size[,2:]
spec_struct.particle_size=particle_size
citation=citation[,2:]
spec_struct.citation=citation
type=type[,2:]
spec_struct.group=type
mineral_number=mineral_number[,2:]
spec_struct.dana_mineral_number=mineral_number
chemical_formula=chemical_formula[,2:]
spec_struct.chemical_formula=chemical_formula
subgroup=subgroup[,2:]
spec_struct.type_subgroup=subgroup
}
###################################
# speclib fields "Spectral Info"(optional)
###################################
if(lib=="yes"){
date=pause("Enter the date that the spectra were acquired e.g. 2006-09-30 \n")
spec_struct.analysis_date= date[1:(strlen(date)-1)]
printf("Enter the spectrometer on which that the spectra were acquired:\n")
printf("1 = Nicolet Nexus 670 FTIR,2= Nicolet Nexus 470 FTIR,\n")
inst=pause("3 = D&P Field Spectrometer, 4 = Mattson-Cygnus 100 FTIR \n")
if(inst[1,1]=="1"){
spec_struct.instrument= "Nicolet Nexus 670 FTIR"
}else{
if(inst[1,1]=="2"){
spec_struct.instrument= "Nicolet Nexus 470 FTIR"
}else{
if(inst[1,1]=="3"){
spec_struct.instrument= "D&P Field Spectrometer Model 100-IR"
}else{
if(inst[1,1]=="4"){
spec_struct.instrument= "Mattson-Cygnus 100 FTIR"
}else{
spec_struct.instrument=inst[1:(strlen(inst)-1)]
}}}}
wn_low=pause("Enter the lower wavenumber range in cm-1\n")
spec_struct.wavenum_range_low = wn_low[1:(strlen(wn_low)-1)]
wn_hi=pause("Enter the upper wavenumber range in cm-1\n")
spec_struct.wavenum_range_high = wn_hi[1:(strlen(wn_hi)-1)]
res=pause("Enter the spectral resolution in cm-1\n")
spec_struct.resolution=res[1:(strlen(res)-1)]
fov=pause("Enter the field of view in centimeters\n")
spec_struct.field_of_view=fov[1:(strlen(fov)-1)]
person=pause("Enter the name of the person who acquired the spectra\n")
spec_struct.spectral_analysis_person=person[1:(strlen(person)-1)]
lab=pause("Enter the source lab for the spectra : 1=ASU, 2=UofH, 3=Brown, 4=Pitt, 5=Indian Inst of Tech, 6=Other \n")
if(lab[1,1]=="1"){
spec_struct.source_lab="Arizona State University"
}else{
if(lab[1,1]=="2"){
spec_struct.source_lab="University of Hawaii"
}else{
if(lab[1,1]=="3"){
spec_struct.source_lab="Brown University"
}else{
if(lab[1,1]=="4"){
spec_struct.source_lab="University of Pittsburg"
}else{
if(lab[1,1]=="5"){
spec_struct.source_lab="Indian Institute of Technology"
}else{
if(lab[1,1]==""){
spec_struct.source_lab="Other"
}}}}}}
}
###################################
# write data to file
###################################
if(deh=="yes") {
printf("\nRunning dehyd2()\n")
spec_struct.orig_data=spec_struct.data
spec_struct.data=dehyd2(spec_struct.orig_data,spec_struct.xaxis)
}
write(spec_struct,directory + out,hdf, force=1)
indiv=struct()
len2=dim(spec_struct.data)[2]
if (rat=="yes"){
for(l=1;l<=len2;l+=1) {
indiv.xaxis=spec_struct.xaxis
indiv.label="emcal emissivity of " + spec_struct.sample_name[,l]
indiv.data=spec_struct.data[,l]
write_vm(indiv,directory + spec_struct.sample_name[,l] + ".txt")
}
}
if (ind=="yes"){
printf("Writing individual files\n")
for(l=1;l<=len2;l+=1) {
indiv.data=spec_struct.data[,l]
indiv.xaxis=spec_struct.xaxis
indiv.label=spec_struct.sample_name[,l]
write_vm(indiv,directory + indiv.label + ".txt")
}
}
return(spec_struct)
}
define dehyd2() {
#added $DV_EX support
#Added $DV_SCRIPT_FILES support - 5-13-08
if($ARGC==0) {
printf("\n\"Dehydrate\" a lab spectrum\n")
printf("A resampled ratio of water sensitive bands is computed and compared to a pure water specturm \n\n")
printf("$1 = emissivity data\n")
printf("$2 = Xaxis\n")
printf("stolen from j.bandfield (c.edwards Oct 10,2006)\n\n")
return(null);
}
dehyd_water_data=$1
xaxis=$2
verbose=0
dehyd_water=read_vm($DV_SCRIPT_FILES+"/dehyd_water.txt")
dehyd_water=resample(dehyd_water.data,dehyd_water.xaxis,xaxis)
verbose=3
wat_index=1-(avg((dehyd_water.data[,,16]//dehyd_water.data[,,26]//dehyd_water.data[,,30]//dehyd_water.data[,,42]//dehyd_water.data[,,43]//dehyd_water.data[,,55]//dehyd_water.data[,,66]//dehyd_water.data[,,68]//dehyd_water.data[,,80])/(dehyd_water.data[,,12]//dehyd_water.data[,,22]//dehyd_water.data[,,32]//dehyd_water.data[,,39]//dehyd_water.data[,,46]//dehyd_water.data[,,57]//dehyd_water.data[,,64]//dehyd_water.data[,,70]//dehyd_water.data[,,78]))+avg((dehyd_water.data[,,634]//dehyd_water.data[,,643]//dehyd_water.data[,,653]//dehyd_water.data[,,679]//dehyd_water.data[,,687]//dehyd_water.data[,,697]//dehyd_water.data[,,706]//dehyd_water.data[,,755]//dehyd_water.data[,,771]//dehyd_water.data[,,780]//dehyd_water.data[,,788]//dehyd_water.data[,,797])/(dehyd_water.data[,,631]//dehyd_water.data[,,645]//dehyd_water.data[,,651]//dehyd_water.data[,,677]//dehyd_water.data[,,683]//dehyd_water.data[,,700]//dehyd_water.data[,,709]//dehyd_water.data[,,758]//dehyd_water.data[,,769]//dehyd_water.data[,,775]//dehyd_water.data[,,785]//dehyd_water.data[,,793])))*0.5
for(i=1;i<=dim(dehyd_water_data)[1];i+=1) {
for(j=1;j<=dim(dehyd_water_data)[2];j+=1) {
index=1-(avg((dehyd_water_data[i,j,16]//dehyd_water_data[i,j,26]//dehyd_water_data[i,j,30]//dehyd_water_data[i,j,42]//dehyd_water_data[i,j,43]//dehyd_water_data[i,j,55]//dehyd_water_data[i,j,66]//dehyd_water_data[i,j,68]//dehyd_water_data[i,j,80])/(dehyd_water_data[i,j,12]//dehyd_water_data[i,j,22]//dehyd_water_data[i,j,32]//dehyd_water_data[i,j,39]//dehyd_water_data[i,j,46]//dehyd_water_data[i,j,57]//dehyd_water_data[i,j,64]//dehyd_water_data[i,j,70]//dehyd_water_data[i,j,78]))+avg((dehyd_water_data[i,j,634]//dehyd_water_data[i,j,643]//dehyd_water_data[i,j,653]//dehyd_water_data[i,j,679]//dehyd_water_data[i,j,687]//dehyd_water_data[i,j,697]//dehyd_water_data[i,j,706]//dehyd_water_data[i,j,755]//dehyd_water_data[i,j,771]//dehyd_water_data[i,j,780]//dehyd_water_data[i,j,788]//dehyd_water_data[i,j,797])/(dehyd_water_data[i,j,631]//dehyd_water_data[i,j,645]//dehyd_water_data[i,j,651]//dehyd_water_data[i,j,677]//dehyd_water_data[i,j,683]//dehyd_water_data[i,j,700]//dehyd_water_data[i,j,709]//dehyd_water_data[i,j,758]//dehyd_water_data[i,j,769]//dehyd_water_data[i,j,775]//dehyd_water_data[i,j,785]//dehyd_water_data[i,j,793])))*0.5
ratio=index/wat_index
fix=1-((1-dehyd_water.data)*ratio)
dehyd_water_data[i,j]=float(dehyd_water_data[i,j]/fix)
}
}
#vm code
# index=1-(avg(dehyd_water_data[,,16//26//30//42//43//55//66//68//80]/dehyd_water_data[,,12//22//32//39//46//57//64//70//78])+avg(dehyd_water_data[,,634//643//653//679//687//697//706//755//771//780//788//797]/dehyd_water_data[,,631//645//651//677//683//700//709//758//769//775//785//793]))*0.5
# wat_index=1-(avg(dehyd_water.data[,,16//26//30//42//43//55//66//68//80]/dehyd_water.data[,,12//22//32//39//46//57//64//70//78])+avg(dehyd_water.data[,,634//643//653//679//687//697//706//755//771//780//788//797]/dehyd_water.data[,,631//645//651//677//683//700//709//758//769//775//785//793]))*0.5
return(dehyd_water_data)
}
define sublib() {
if($ARGC==0) {
printf("Subset a spectral library structure\n\n")
printf("$1 = library structure\n")
printf("$2 = subet() string\n\n")
printf("NOTE: This function uses subset(). See the subset() help for more information\n\n")
printf("C.Edawrds 4-8-11\n\n")
return(null)
}
lib=$1
if($ARGC<2) {
subsetstr=""
} else {
subsetstr=$2
}
sublib={}
key=get_struct_key(lib)
for(i=1;i<=length(lib);i+=1) {
data=lib[i]
if(type(data)!="TEXT") {
if(dim(data)[1]==1 && dim(data)[2]!=1) {
tmp=subset(data,y=strsub(subsetstr,",",""))
printf("Data in X\n")
} else if (dim(data)[1]!=1 && dim(data)[2]==1) {
tmp=subset(data,x=strsub(subsetstr,",",""))
printf("Data in y\n")
} else if (dim(data)[1]!=1 && dim(data)[2]!=1) {
tmp=subset(data,subsetstr)
} else {
tmp=data
}
} else {
tmp=subset(data,y=strsub(subsetstr,",",""))
printf("Data is TEXT\n")
}
add_struct(obj=sublib,name=key[,i],value=tmp)
}
return(sublib)
}
define sma(wave1,wave2,surface,bb,band,group,forceall,exclude,sort,noerr,resample,notchco2,nn) {
if ($ARGC < 2 || $ARGC > 3) {
printf (" \n")
printf (" sma=spectral mixture analysis \n")
printf (" \n")
printf (" $1 = mixed cube (or structure with .data and .xaxis elements)\n")
printf (" $2 = endmembers. $2 can be a standard structure with at least the following \n")
printf (" elements: .data(arranged [#specs-x, #specs-y, #bands]), .label, \n ")
printf (" and .xaxis (arranged [1,1,x-axis values]) or may be an array (arranged \n ")
printf (" [#specs-x, #specs-y, #bands] \n")
printf (" \n")
printf (" $1 and $2 must have an equal number of spectral bands. \n")
printf (" \n")
printf (" [$3] = endmembers to be forced despite negative concentrations (optional). The \n")
printf (" forced library must be in the same format as $2 \n")
printf (" elements. \n")
printf (" \n")
printf ("Options: \n")
printf (" surface =1 for atm-removed and modeled surface spectra in addition to normal output \n")
printf (" wave1 = wavenumber/micron/channel value for modeling (low); default=400 cm-1 \n")
printf (" wave2 =wavebumber/micron/channel value for modeling (high); default=1600 cm-1 \n")
printf (" band = 1 if wave1 and wave2 are given in channel numbers, instead of frequency (cm-1) \n")
printf (" or wavelength \n")
printf (" bb = 0 do not add blackbody to endmember matrix; bb added by default \n")
printf (" group = 1 additional output of model-derived modes summed by mineral group (eg. feldspar \n ")
printf (" abundance, pyroxene abundance). $2 must contain .group. Default is 1\n ")
printf (" \n")
printf (" forceall = 1 to keep negative endmember percentages (requires nn = 0 as well)\n")
printf (" exclude= vector list of endmembers to exclude (ex:2//39//40). \n")
printf (" noerr=1 do not calculate statistical errors (recommended for large cubes) \n")
printf (" resample=1 resample the mixed cube structure to the endmember library (default=0)\n")
printf (" resample=2 resample the endmember library to the mixed cube structure (default=0)\n")
printf (" sort = 1 if you want abundances sorted for each pixel; default is sorting \n")
printf (" notchco2 = 1 automatically remove CO2 band in data/library/forced centered at 669 +/- 74 wn\n")
printf (" notchco2 = N removve CO2 band in data/library/forced centered at 669 +/- N wn (Default is 0)\n")
printf (" nn = 0 (default) to use old fitting algorithm\n")
printf (" nn = 1 to use algorithm optimized for nonnegative least squares\n")
printf (" (NNLS, outlined in Lawson & Hanson 1974), which is faster and more accurate\n")
printf ("\n\n")
printf (" Examples: \n")
printf (" minlib=read('my_mineral_lib.hdf') \n")
printf (" atmlib=read('my_atm_lib.hdf') \n")
printf (" nw_scene=read('mtes_nw_raster.hdf') \n")
printf (" \n")
printf (" final=sma(nw_scene.spec, minlib, atmlib,wave1=230, wave2=1400,surface=1) \n")
printf (" --> Unmix the entire scene with mineral and atm endmembers, between 230-1400 cm-1, \n")
printf (" return normalized concentrations and atm-removed emissivity in addition to normal output. \n")
printf (" \n")
printf (" final2=sma(nw_scene.spec[1:30,1:30,], minlib, band=1, wave1=3, wave2=160) \n")
printf (" --> Unmix the top 30x30 spectra in the cube with mineral endmembers only, between channels \n")
printf (" 3-160. \n")
printf (" \n")
printf (" See also sort_cube, sum_group_conc, summary_sma, and summary_sma_group \n")
printf (" \n")
return(0)
}
# deanne rogers 12/2007
# modified to handle resample option 5/16/08 cedwards
# modified to add notchco2 option from deconvolve 7/30/10
# Modified by S. Marshall, 11-23-2010 (added nn)
# Modified to check for spectral libraries in the correct direction (for some reason that is x as it completely fails if its in y) 2-3-2011
# Read in emissivity spectra and library(ies)
em_cube=$1
libstruc=$2
if(HasValue(sort)==0) {
sort=1
}
if(HasValue(surface)==0) {
surface=0
}
#check for groups
if(type(libstruc)=="STRUCT") {
if(HasValue(get_struct(libstruc,"group"))) {
if(HasValue(group)==0) {
group=1
}
} else {
if(HasValue(group)==0) {
group=1
}
}
} else {
if(HasValue(group)==0) {
group=0
}
}
# Create output structure
smaout=struct()
#The spectral library ouputs spectra in a fashion that is not useable by sma (so make it useable)
if(type(libstruc)=="STRUCT") {
if(HasValue(get_struct(libstruc,"sample_id")) && HasValue(get_struct(libstruc,"shortlabel"))==0){
libstruc=reform_speclib(libstruc)
}
#unfortunately now we check to make sure the data is in the x direction for the library.
if(dim(libstruc.data)[1]==1 && dim(libstruc.data)[2]>=1) {
libstruc.data=translate(libstruc.data,y,x)
}
} else {
if(dim(libstruc)[1]==1 && dim(libstruc)[2]>=1) {
libstruc=translate(libstruc,y,x)
}
}
if($ARGC==3) {
atmstruc=$3
if(type(atmstruc)=="STRUCT") {
if(HasValue(get_struct(atmstruc,"sample_id")) && HasValue(get_struct(atmstruc,"shortlabel")) == 0) {
atmstruc=reform_speclib(atmstruc)
}
#unfortunately now we check to make sure the data is in the x direction for the library.
if(dim(atmstruc.data)[1]==1 && dim(atmstruc.data)[2]>=1) {
atmstruc.data=translate(atmstruc.data,y,x)
}
} else {
if(dim(atmstruc)[1]==1 && dim(atmstruc)[2]>=1) {
atmstruc=translate(atmstruc,y,x)
}
}
}
#check to see if we can resample
verbose=0
if(HasValue(resample)==0) resample=0
if((resample==1 || resample==2) && (type(em_cube)!="STRUCT" || HasValue(get_struct(em_cube,"data"))==0 || HasValue(get_struct(em_cube,"xaxis"))==0)) {
printf("A mixed structure with .data and .xaxis elements is requred for the \"resample\" option\n")
return(null)
}
#if the em_cube is a struct the we can possibly resample
if(type(em_cube)=="STRUCT") {
#set up labels for summary_sma(_group)
if(HasValue(get_struct(em_cube,"sample_name"))) smaout.sample_name=em_cube.sample_name
if(HasValue(get_struct(em_cube,"label"))) smaout.label=em_cube.label
#resample the spectra to the endmember library
if(resample==1 && HasValue(get_struct(em_cube,"xaxis")) && HasValue(get_struct(em_cube,"data")) && HasValue(get_struct(libstruc,"xaxis"))) {
smaout.original={}
em_cube.data=float(em_cube.data)
em_cube.xaxis=float(em_cube.xaxis)
smaout.original.data=em_cube.data
smaout.original.xaxis=em_cube.xaxis
if($DV_OS=="win") {
em_cube.data=resampleudf(em_cube.data,em_cube.xaxis,libstruc.xaxis).data
} else {
em_cube.data=resample(em_cube.data,em_cube.xaxis,libstruc.xaxis).data
}
} else if (resample==1 && (HasValue(get_struct(em_cube,"xaxis"))==0 || HasValue(get_struct(em_cube,"data"))==0 || HasValue(get_struct(libstruc,"xaxis"))==0)){
printf("Mixed cube strucure with .xaxis and .data elements is required for the \"resample=1\" option\n")
printf("Additionally an endmember strucutre with a .xaxis element is required\n")
verbose=1
return(null)
}
#resample the endmember library to the spectra
if(resample==2 && HasValue(get_struct(em_cube,"data")) && HasValue(get_struct(libstruc,"data")) && HasValue(get_struct(libstruc,"xaxis"))) {
libstruc.original={}
em_cube.data=float(em_cube.data)
em_cube.xaxis=float(em_cube.xaxis)
libstruc.original.data=libstruc.data
libstruc.original.xaxis=libstruc.xaxis
if($DV_OS=="win") {
libstruc.data=resampleudf(libstruc.data,libstruc.xaxis,em_cube.xaxis).data
} else {
libstruc.data=resample(libstruc.data,libstruc.xaxis,em_cube.xaxis).data
}
libstruc.xaxis=em_cube.xaxis
} else if(resample==2 && (HasValue(get_struct(em_cube,"data"))==0 || HasValue(get_struct(libstruc,"data"))==0 || HasValue(get_struct(libstruc,"xaxis"))==0)) {
printf("Mixed cube strucure with .xaxis and .data elements is required for the \"resample=2\" option\n")
printf("Additionally an endmember structure with .xaxis and .data elements is required\n")
verbose=1
return(null)
}
if(HasValue(get_struct(em_cube,"data"))){
em_cube=em_cube.data
} else if (HasValue(get_struct(em_cube,"data"))==0){
printf("No .data element was found in the mixed structure. Unable to continue.\n")
verbose=1
return(null)
}
} else {
em_cube=float(em_cube)
}
verbose=1
# Set flags for later use
labelflag=1
grpflag=1
# Determine if library is standard structure
if (type(libstruc)=="STRUCT"){
libstruc.data=float(libstruc.data)
if (HasValue(libstruc.data)==1) {
endmem=libstruc.data
libchan=dim(endmem)
} else {
endmem=libstruc
libchan=dim(endmem)
}
# Check if library has labels
if (HasValue(libstruc.label)==1) {
names=libstruc.label
} else{
names=0
labelflag=0
}
# Check if library has xaxis
if (HasValue(libstruc.xaxis)==1) {
x_ax=libstruc.xaxis
} else {
x_ax=create(1,1,libchan[3],start=1,step=1,format=float,org=bsq)
}
# If group=1 option was selected, check if library has groups
if(group==1) {
if (HasValue(libstruc.group)==1){
grpnames=libstruc.group
} else {
grpnames=0
grpflag=0
}
}
} else {
libstruc=float(libstruc)
endmem=libstruc
libchan=dim(endmem)
x_ax=create(1,1,libchan[3],start=1,step=1,format=float,org=bsq)
names=0
labelflag=0
grpnames=0
grpflag=0
if (HasValue(band)==0) {
wave1=1
wave2=dim(em_cube)[3]
printf("No wavenumber x-axis found. Using full spectral range. \n")
}
}
# Check that library and em_cube have the same number of channels
dimend=dim(em_cube)
if (dimend[3] != libchan[3]) {
printf(" $1 and $2 must have the same z-dimension. \n")
return(0)
}
if(sum(sum(endmem,axis=z)==0)!=0) {
printf("\nOne of the endmembers in your endmember library is composed of all zeros\n")
printf("Please check your library and try again\n\n")
return(null)
}
#Here is an auto/custom crop for the co2 band
if(HasValue(notchco2)==0){
notchco2=0
} else {
co2wn=669
if(notchco2==1) {
# use default value of 7 (historical)
# was 7 tes bands and is now stored as equivalent mini-TES wavenumbers
wco2 = 74.06
} else {
wco2 = notchco2
}
co2chan=int(round(interp(create(1,1,dim(x_ax)[3],start=1),x_ax,co2wn)))
#scale wavenumbers to channels
dwavex=x_ax[,,co2chan+1]-x_ax[,,co2chan]
nco2 = int(round(wco2/dwavex))
co2low = co2chan - nco2
co2high = co2chan + nco2
# emissivity cube
etmp = cat(em_cube[,,1:co2low], em_cube[,,co2high:], z)
em_cube = etmp
# library spectra
mtmp = cat(endmem[,,1:co2low], endmem[,,co2high:], z)
endmem = mtmp
# xaxis
xtmp = cat(x_ax[,,1:co2low], x_ax[,, co2high:], z)
x_ax = xtmp
if(type(libstruc)=="STRUCT") {
if(HasValue(libstruc.xaxis)==1) {
libstruc.xaxis=x_ax
libstruc.data=endmem
}
}
if(type(em_cube)=="STRUCT") {
if(HasValue(em_cube.xaxis)) {
em_cube.xaxis=x_ax
em_cube.data=em_cube
}
}
smaout.notchco2={}
smaout.notchco2.wco2=wco2
smaout.notchco2.co2chanlow=co2low
smaout.notchco2.co2chanhigh=co2high
smaout.notchco2.dwavex=dwavex
smaout.notchco2.co2low=co2wn-wco2
smaout.notchco2.co2high=co2wn+wco2
libchan=dim(endmem)
dimend=dim(em_cube)
}
#setup the output xaxis for easy plotting later.
smaout.xaxis=x_ax
#Check if there is an atmospheric library
if ($ARGC==3) {
smaout.forcedlib=atmstruc
#Check if atmospheric library is standard structure
if(type(atmstruc)=="STRUCT"){
if (HasValue(atmstruc.data)==1) {
norem=atmstruc.data
} else {
norem=atmstruc
}
#Check that the atm library has an xaxis
if (HasValue(atmstruc.xaxis)==1) {
atmxax=atmstruc.xaxis
} else {
atmxax=create(1,1,dim(norem)[3],start=1,step=1,format=float,org=bsq)
}
#Check that the atm library has labels
if (HasValue(atmstruc.label)==1) {
atmnames=atmstruc.label
} else {
atmnames=0
if (labelflag!=0) {
printf(" $2 and $3 must both have labels or both not have labels. \n")
}
}
# Check that the atm library has group names
if(group==1) {
if (HasValue(atmstruc.group)==1) {
atmgrpnames=atmstruc.group
} else {
atmgrpnames=0
if (grpflag!=0) {
printf(" $2 and $3 must both have group labels or both not have group labels. \n")
return(0)
}
}
}
} else {
#if the atmstruc is not a structure, we go by index only
norem=atmstruc
atmxax=create(1,1,dim(norem)[3],start=1,step=1,format=float,org=bsq)
atmnames=0
if (labelflag!=0) {
printf(" $2 and $3 must both have labels or both not have labels. \n")
}
atmgrpnames=0
if (grpflag!=0) {
printf(" $2 and $3 must both have group labels or both not have group labels. \n")
}
}
#again with the auto/custom co2 band removal (this time for the forcedlib)
if(notchco2!=0) {
co2chan=int(round(interp(create(1,1,dim(atmxax)[3],start=1),atmxax,co2wn)))
#scale wavenumbers to channels
dwavex=atmxax[,,co2chan+1]-atmxax[,,co2chan]
nco2 = int(round(wco2/dwavex))
co2low = co2chan - nco2
co2high = co2chan + nco2
# emissivity cube
noremtmp = cat(norem[,,1:co2low], norem[,,co2high:], z)
norem=noremtmp
# xaxis
atmxaxtmp = cat(atmxax[,,1:co2low], atmxax[,, co2high:], z)
atmxax=atmxaxtmp
if(type(smaout.forcedlib)=="STRUCT") {
if(HasValue(smaout.forcedlib.xaxis)==1) {
smaout.forcedlib.xaxis=atmxax
smaout.forcedlib.data=norem
}
}
}
#Check that the two libraries have the same number of channels
if (libchan[3] != dim(norem)[3]) {
printf(" $2 and $3 must have the same z-dimension. \n")
return(0)
}
#Cat atmospheric endmembers onto libst
endmem=cat(endmem,norem,x)
#Determine number of forced lib spectra
forced_specs=dim(norem)[1]
} else {
forced_specs=0
smaout.forcedlib=struct()
smaout.forcedlib.label="none"
}
# Remember full size of library so that it can be expanded back to original size later
dimfull=dim(endmem)
# Reduce library if exclude option is set (set excluded endmembers to zero)
if (HasValue(exclude)) {
smaout.excluded=int(translate(unique(sort(exclude)),y,x))
temp_end=endmem[1]
countv=1
count=1
while(count <= dim(endmem)[1]) {
while(countv<=dim(smaout.excluded)[1]) {
if(count==smaout.excluded[countv]) {
countv=countv+1
} else {
temp_end=cat(temp_end,endmem[count],x)
}
count=count+1
}
if(count <= dim(endmem)[1]) {
temp_end=cat(temp_end,endmem[count],x)
}
count=count+1
}
endmem=temp_end[2:]
}
# Determine number of endmembers not to be removed, add blackbody spectrum to library if requested
# There was previously a problem with this part
# Flat blackbody spectrum was added to endmem if bb omitted, but not if bb=1
# Fixed 11/23/2010
if (HasValue(bb)==1) {
if (equals(bb, 0) || equals(bb, 1)) { # Now checking to ensure that bb is either 0 or 1
bbval=bb
} else {
printf("Error! Value for bb must be 0 or 1!\n")
return(byte(0))
}
} else {
bbval=1
#endmem=cat(endmem,create(1,1,dim(endmem)[3],format=float,start=1,step=0),x) # This belongs outside
}
if (bbval == 1) {
endmem=cat(endmem,clone(float(1), 1,1,dim(endmem)[3]),x)
}
num_no_rem=forced_specs+bbval
# Check if errors should be calculated
if(HasValue(noerr)==1) {
no_calc_error=1
}else{
no_calc_error=0
}
# Determine wavelength range from user input
if (HasValue(wave1)) {
wv1=wave1
} else {
wv1=400
}
if (HasValue(wave2)){
wv2=wave2
} else {
wv2=1600
}
#Crop xaxis according to user input
if (HasValue(band)==1){ # if wave1 and wave 2 are input as channels
stband=wv1
endband=wv2
new_xax=x_ax[,,stband:endband]
new_endmem=endmem[,,stband:endband]
new_cube=em_cube[,,stband:endband]
} else {
# Convert wavenumber or micron start/stop input into channel/index start/stop bands
if (wv1 < x_ax[,,1]) {
stband=1
} else {
count=1
while(x_ax[,,count] < wv1) {
count=count+1
}
stband=count
}
if (x_ax[,,libchan[3]] <= wv2) {
endband=libchan[3]
} else {
count=1
while(x_ax[,,libchan[3]-count] > wv2) {
count=count+1
}
endband=libchan[3]-count
}
if (endband < stband) { #If the library xaxis (e.g., THEMIS) is backwards, switch the start and end channels
tmpband=stband
stband=endband
endband=tmpband
}
new_xax=x_ax[,,stband:endband]
new_endmem=endmem[,,stband:endband]
new_cube=em_cube[,,stband:endband]
}
# Determine whether to use optimimal nonnegative least squares algorithm or original sma algorithm
# Added 11/23/2010
if (HasValue(nn) == 1) {
if (equals(nn, 0) || equals(nn, 1)) { # Now checking to ensure that nn is either 0 or 1
nnval = nn
} else {
printf("\nError! Value for nn must be 0 or 1!\n\n")
return(byte(0))
}
} else {
nnval = 0 # By default, it will NOT use the nnls algorithm
}
# Transpose endmember matrix for matrix operations
new_endmem=translate(new_endmem,y,z)
new_xax=translate(new_xax,y,z)
# Get dimensions of emissivity cube and endmembers
mixvecs=dim(new_cube)
cnt=dim(new_endmem)
nend=cnt[1]
nsamp=(cnt[2])
# Warn user of underdetermined solution
if (nend > nsamp) {
printf ("You have more endmembers (%i) than bands (%i). \n",nend,nsamp)
return(0)
}
# Take strips off the cube and stack in the y-direction
em=new_cube[1,,]
em=translate(em, from=z, to=x)
count=2
for (count=2; count<=mixvecs[1,,]; count+=1){
# while (count <= mixvecs[1,,])
# em_temp=translate(new_cube[(count),,],from=x, to=y)
em_temp=translate(new_cube[(count),,],from=z, to=x)
em=cat(em,em_temp,y)
# count=count+1
}
emd = dim(em)
#create empty error matrix for populating
st_err=clone(float(0), nend,emd[2])
# Finally, to the core. Create transpose and square endmember matrix.
endmemT=translate(new_endmem,x,y)
emT=translate(em, from=x, to=y)
square=mxm(endmemT,new_endmem)
covm = byte(0) # Giving it an initial value, to avoid errors later if no_calc_error != 0
# Core section (different depending on value of nn) begins here
if (nnval == 0) { # Will use original algorithm
# Invert endmember matrix
square=minvert(square)
# Multiply inverse square matrix by endmember matrix.
sq_end=mxm(new_endmem,square)
# Multiply sq_end by array of emiss spectra producing the percentage matrix
perc=mxm(em,sq_end)
# Detect mixed spectra with neg. endmember percentages, loop through each spectrum
npositive=dim(new_endmem)
npositive=(npositive[1]-num_no_rem)
} else { # If nn == 1
# Main block
cnt1 = 1 # Counter to iterate through each (mixed) test spectrum
perc = clone(double(-32768), emd[2], nend)
if (num_no_rem == 0) { # If there's no blackbody component to be included, this is easy
for (cnt1 = 1; cnt1 <= emd[2]; cnt1++) {
perc[cnt1,] = lsqnn(new_endmem, emT[cnt1,], 0)
}
} else { # Need to be careful with forced endmembers
# Must call lsqsn to handle forced endmembers
for (cnt1 = 1; cnt1 <= emd[2]; cnt1++) {
perc[cnt1,] = lsqsn(new_endmem, emT[cnt1,], num_no_rem, 0)
}
}
perc=translate(perc,x,y) # So that it's in the proper shape for the next loop
}
if (HasValue(forceall)==0 || nnval!=0) { #If you don't want negative percentages, do this loop. Otherwise, skip to the end.
# If forceall == 1 and nn == 1, it has already used the nnls algorithm, and all concentrations will be positive
cnt1=1
while (cnt1 <= emd[2]) {
if (nnval == 0) {
# Compose library w/ no negative endmembers. Loop until no negative surface concentrations.
iterflag=0
while (min(perc[1:npositive,cnt1])<0) {
# Find positive percentages and add to new endmember matrix
newlib=new_endmem[1]
cnt2=1
while (cnt2 <= npositive) {
if (perc[cnt2,cnt1] > 0) {
newlib=newlib//new_endmem[cnt2]
}
cnt2=cnt2+1
}
if (bbval==1 || num_no_rem>0) {
newlib=newlib//new_endmem[(npositive+1):(npositive+num_no_rem)]
}
newlib=newlib[2:]
# Create transpose and square new endmember matrix
newlibT=translate(newlib,x,y)
square=mxm(newlibT,newlib)
# Invert new endmember matrix
square=minvert(square)
# Multiply inverse square matrix by new endmember matrix
sq_end=mxm(newlib,square)
# Multiply sq_end by transpose of emiss spectra producing the new 1x?x1 perc matrix
newperc=mxm(em[,cnt1],sq_end)
# Set negative percentages to zero before replacing new iteration percentages
negcheck=perc[1:(npositive),cnt1]
negcheck[where(negcheck<0)]=0
perc[1:(npositive),cnt1]=negcheck
# Replace original percentages with new percentages
cnt5=1
cnt6=1
while (cnt5 <= npositive) {
if (perc[cnt5,cnt1] > 0) {
perc[cnt5,cnt1]=newperc[cnt6]
cnt6=cnt6+1
}
cnt5=cnt5+1
}
#Add blackbody percentage back to new percent matrix
if (bbval==1 || num_no_rem > 0) {
perc[(npositive+1):(npositive+num_no_rem),cnt1]=newperc[cnt6:]
}
iterflag=1
}
#If no iterations were needed, set newlib and newperc to new_endmem and perc[,cnt1]
if(iterflag==0) {
newlib=new_endmem
newperc=perc[,cnt1]
}
} else { # If nn == 1
# Filling in needed variables
# Must avoid instances where A[i] == 0
#The past definition of npositive was incorrect.
#It is just the number of endmembers minus the number of forced endmembers
#It is not the number of positive endmemebers...Its the number that can only be positive or zero....
npositive = int(dim(perc[1:(nend-num_no_rem),cnt1])[1])
newlib = new_endmem[1,,]*.0 # Have to give it an initial value
newperc = perc[1,cnt1,]*0.
cnt2 = 1 # Counter
for (cnt2 = 1; cnt2 <= (nend - num_no_rem); cnt2++) {
if (perc[cnt2,cnt1,] > 0.) {
newlib = newlib//new_endmem[cnt2,,]
newperc = newperc//perc[cnt2,cnt1,]
}
}
if(dim(newlib)[1]!=1) {
newlib = newlib[2:,,] # Remove initial column
newperc = newperc[2:,,]
}
if (bbval==1 || num_no_rem>0) { # I think this is redundant, since num_no_rem includes the blackbody
newlib=newlib//new_endmem[(nend-num_no_rem+1):nend,,]
newperc = newperc//perc[(nend-num_no_rem+1):nend,cnt1,]
}
}
#Generate covariance matrix
if(no_calc_error==0) {
emT=translate(em[,cnt1,1],x,y)
newpercT=translate(newperc,x,y)
newlibT=translate(newlib,x,y)
term1=mxm(em[,cnt1,1],emT)
term2=mxm(newperc,newlibT)
term3=mxm(term2,emT)
dimnewlib=dim(newlib)
term4=nsamp-dimnewlib[1]
term5=minvert(mxm(newlibT,newlib))
covm=((term1-term3)/term4)*(term5)
#take statistical errors from sqrt of diagonal of covariance matrix
tmp_err=create(1,1,dimnewlib[1],format=float,start=0,step=0)
for (count=1; count<=dim(newlib)[1];count+=1) {
tmp_err[,,count]=sqrt(abs(covm[count,count]))
}
#expand st_err into same length as perc
count=1
for (i=1; i<=dim(perc)[1]; i+=1) {
if (perc[i,cnt1] != 0) {
st_err[i,cnt1]=tmp_err[,,count]
count++
}
}
}
#end error calculations
cnt1++
}
}
# End of core section (should be same after this for either value of nn)
smaout.covm=covm
# Reconstruct original spectra
perc=translate(perc,x,y)
endmem=translate(endmem,y,z)
recon=mxm(new_endmem,perc) #For RMS error calculation
reconfull=mxm(endmem,perc) #For modeled spectrum
dimperc=dim(perc)
# Reconstruct modeled surface spectra
if (surface==1) {
if (bbval==1) {
endmem_modsur=endmem[1:(npositive),,]//endmem[(nend),,]
perc_modsur=cat(perc[,1:(npositive),],perc[,(nend),],axis=y)
}
if (bbval==0){
endmem_modsur=endmem[1:(npositive),,]
perc_modsur=perc[,1:(npositive),]
}
modsur=mxm(endmem_modsur,perc_modsur)
}
# Reconstruct remaining data cubes from output
recon_cube=create(x=dimend[1,,],y=dimend[2,,],z=dimend[3,,],format=float,start=0,step=0)
perc_cube=create(x=mixvecs[1,,],y=mixvecs[2,,],z=dimperc[2,,], format=float,start=0,step=0)
if(no_calc_error==0) {
err_cube=create(x=mixvecs[1,,],y=mixvecs[2,,],z=dimperc[2,,], format=float,start=0,step=0)
}
#If your input mixed cube has a y-dimension > 1
if (mixvecs[2] > 1) {
perc=translate(perc,x,y)
recon=translate(recon,x,y)
reconfull=translate(reconfull,x,y)
for (count=1; count<=mixvecs[1,,]; count+=1){
start_spec=(count-1)*mixvecs[2,,]+1
stop_spec=count*mixvecs[2,,]
recon_temp=translate(reconfull[,start_spec:stop_spec,],from=x,to=z)
recon_cube[count,,]=recon_temp
perc_temp=translate(perc[,start_spec:stop_spec,],from=x,to=z)
perc_cube[count,,]=perc_temp
if(no_calc_error==0){
err_temp=translate(st_err[,start_spec:stop_spec,],from=x,to=z)
err_cube[count,,]=err_temp
}
}
} else {
recon_cube=translate(reconfull, from=y, to=z)
perc_cube=translate(perc, from=y, to=z)
if(no_calc_error==0) {
err_cube=translate(st_err,from=y,to=z)
err_cube=translate(err_cube,from=x,to=z)
}
}
# Generate blackbody plane
# There was a problem with this section where bb=0; fixed 11/23/2010
if (bbval == 1) {
smaout.bb = float(perc_cube[,,nend])
if (no_calc_error == 0) {
smaout.bberror = float(err_cube[,,nend])
err_cube = err_cube[,,1:(nend-1)]
}
perc_cube = perc_cube[,,1:(nend-1)]
smaout.normconc = float(perc_cube/sum(perc_cube, axis=z)*100.0)
} else { # Need this for sum_group_conc
smaout.bb = float(0*perc_cube[,,1])
if (no_calc_error == 0) {
smaout.bberror = smaout.bb
#err_cube=err_cube[,,1:(nend-1)]
}
}
#Calculate RMS error
rms_plane=new_cube-recon_cube[,,stband:endband]
rms_planefull=em_cube-recon_cube
# Create rematm and modsur cubes if requested
if(surface==1) {
modsur_cube=create(x=dimend[1],y=dimend[2],z=dimend[3],format=float,start=0,step=0)
modsur=translate(modsur,x,y)
for(count=1; count<=mixvecs[1,,]; count+=1){
start_spec=(count-1)*mixvecs[2,,]+1
stop_spec=count*mixvecs[2,,]
modsur_temp=translate(modsur[,start_spec:stop_spec,],from=x,to=z)
modsur_cube[count,,]=modsur_temp
}
rematm_cube=modsur_cube+rms_planefull
smaout.modsur=float(modsur_cube)
smaout.rematm=float(rematm_cube)
# Add constant to surface spectra to account for atm percentages
if(bbval==1) {
rematm_cube=rematm_cube+sum(perc_cube[,,(npositive+1):(npositive+num_no_rem-1)],axis=z)
modsur_cube=modsur_cube+sum(perc_cube[,,(npositive+1):(npositive+num_no_rem-1)],axis=z)
smaout.rematm=float(rematm_cube)
smaout.modsur=float(modsur_cube)
}
if(bbval==0) {
rematm_cube=rematm_cube+sum(perc_cube[,,(npositive+1):(npositive+num_no_rem)],axis=z)
modsur_cube=modsur_cube+sum(perc_cube[,,(npositive+1):(npositive+num_no_rem)],axis=z)
smaout.rematm=float(rematm_cube)
smaout.modsur=float(modsur_cube)
}
}
#Back to calculating RMS error
rms_plane=pow(rms_plane,2)
rms_plane=sum(rms_plane,z)
rms_plane=rms_plane/nsamp
rms_plane=pow(rms_plane, 0.5)
# Sort concentrations for each pixel if requested (for plot_sma)
if(sort==1) {
sorted_cube=sort_cube(float(perc_cube))
smaout.sort=int(sorted_cube)
}
# Expand percentage and error cubes to be original size of library (if vectors were excluded)
if(HasValue(exclude)==1) {
tmp_normperc_cube=tmp_perc_cube=create(x=mixvecs[1],y=mixvecs[2],z=dimfull[1],format=float,org=bsq,start=0,step=0)
if(no_calc_error==0){
tmp_err_cube=create(x=mixvecs[1],y=mixvecs[2],z=dimfull[1],format=float,org=bsq,start=0,step=0)
}
#expand out the sort array though...this is tricky
if(sort==1) {
tmp_sort_cube=create(x=mixvecs[1],y=mixvecs[2],z=dimfull[1],format=int,org=bsq,start=0,step=0)
tmp_sort_cube[,,:dim(smaout.sort)[3]]=smaout.sort
for(i=1;i<=dim(smaout.excluded)[1];i+=1) {
tmp_sort_cube[where tmp_sort_cube>=smaout.excluded[i] && tmp_sort_cube!=0]=tmp_sort_cube+1
}
}
count=1
countt=1
counttt=1
while(counttt<=dim(perc_cube)[3]) {
while(countt<=dim(smaout.excluded)[1]) {
if(count==smaout.excluded[countt]) {
countt=countt+1
count=count+1
} else {
tmp_perc_cube[,,count]=perc_cube[,,counttt]
tmp_normperc_cube[,,count]=smaout.normconc[,,counttt]
if(no_calc_error==0) {
tmp_err_cube[,,count]=err_cube[,,counttt]
}
count=count+1
counttt=counttt+1
}
}
if(counttt<=dim(perc_cube)[3]) {
tmp_perc_cube[,,count]=perc_cube[,,counttt]
tmp_normperc_cube[,,count]=smaout.normconc[,,counttt]
if(no_calc_error==0) {
tmp_err_cube[,,count]=err_cube[,,counttt]
}
}
count=count+1
counttt=counttt+1
}
perc_cube=tmp_perc_cube
smaout.normconc=tmp_normperc_cube
if(no_calc_error==0) {
err_cube=tmp_err_cube
}
if(sort==1) {
smaout.sort_full=tmp_sort_cube
}
} else {
if(sort==1) {
smaout.sort_full=smaout.sort
}
}
# Add to the output structure
smaout.measured=float(em_cube)
smaout.modeled=float(recon_cube)
smaout.conc=float(perc_cube)
smaout.rms=float(rms_plane)
smaout.nsamples=nsamp
smaout.endlib=libstruc
if(no_calc_error==0) {
smaout.error=float(err_cube)
smaout.normerror=float(smaout.error[]/sum(smaout.conc[],axis=z))*100
}
# if the input library was not a structure, add .label and .group with "none"
if (labelflag==0) {
smaout.endlib=struct()
smaout.endlib.data=libstruc
smaout.endlib.label="none"
smaout.endlib.group="none"
}
#setup the spectral range ouput
smaout.spectral_range = sprintf("Channels %i-%i", stband,endband)
smaout.spectral_channel_low=stband
smaout.spectral_channel_high=endband
smaout.wave1 = wv1
smaout.wave2 = wv2
#fix for weirdly high (but scaled properly) atm and modled surface spectra
if(surface==1) {
corr=1-(smaout.bb+sum(smaout.conc,axis=z))
smaout.modsur=smaout.modsur+corr
smaout.rematm=smaout.rematm+corr
}
# Group concentrations if requested--set up structure for call to sum_groups
if (group==1) {
if(grpflag==0){
printf("The library does not have .group. Concentrations will not be grouped. \n")
} else {
if($ARGC==3){
smaout.grouped=sum_group_conc(smaout,smaout.forcedlib)
} else {
smaout.grouped=sum_group_conc(smaout)
}
#find the posisition of the "Black body" element
pos=maxchan(smaout.grouped.grouped_labels=="Black body")
#extract the black body
smaout.grouped.grouped_bb=smaout.grouped.grouped_conc[,,pos]
#whittle down the error array and extract the black body error
if(no_calc_error==0) {
smaout.grouped.grouped_bberror=smaout.grouped.grouped_error[,,pos]
smaout.grouped.grouped_error=smaout.grouped.grouped_error[,,:pos-1]
}
#whittle down the conc and labels array
smaout.grouped.grouped_conc=smaout.grouped.grouped_conc[,,:pos-1]
smaout.grouped.grouped_labels=smaout.grouped.grouped_labels[,:pos-1]
#create a grouped_sort for the concentrations
smaout.grouped.grouped_sort_full=sort_cube(smaout.grouped.grouped_conc)
#create a normalized concentration
smaout.grouped.grouped_normconc=float(smaout.grouped.grouped_conc*0)
newsum=sum(smaout.grouped.grouped_conc,axis=z)
smaout.grouped.grouped_normconc=float(smaout.grouped.grouped_conc/newsum*100)
if(no_calc_error==0) {
#create a normalized error array
smaout.grouped.grouped_normerror=float(smaout.grouped.grouped_error[]/sum(smaout.grouped.grouped_conc[],axis=z))*100
}
}
}
if (nnval == 0) {
smaout.algorithm = "Original least squares"
} else {
smaout.algorithm = "Nonnegative least squares"
}
return(smaout)
}
define sort_cube() {
#sort in the z-direction, on a pixel-by-pixel basis
if ($ARGC ==0) {
printf ("Sort a cube in the z-direction, on a pixel-by-pixel basis. \n")
printf (" $1 = input cube \n")
return(0)
}
cube=$1
xdim=dim(cube)
ydim=dim(translate(cube,y,x))
zdim=dim(translate(cube,z,x))
band_list=create(1,1,zdim[1],format=format(cube), start=1,step=1)
sortband_cube=create(xdim,ydim,zdim,format=format(cube), start=0,step=0)
x_cnt=1
y_cnt=1
while (x_cnt <= xdim[1,,]) {
while (y_cnt <= ydim[1,,]) {
strip=create(1,1, zdim[1], format=format(cube))
strip[1,1,]=cube[x_cnt,y_cnt,]
sortband=create(1,1,zdim[1],format=format(cube), start=0,step=0)
two_col=cat(strip, band_list,x)
strip=sort(strip)
amax=max(strip)
count=1
cnt2=1
while (cnt2 <= zdim[1,,]) {
if (count > zdim[1,,]) {
count=1
}
while (count <= zdim[1,,]) {
if (two_col[1,1,count] == amax) {
sortband[1,1,cnt2]=two_col[2,1,count]
if (amax == 0) {
sortband[1,1,cnt2]=0
}
cnt3=zdim[1]-cnt2
cnt2=cnt2+1
strip=strip[1,1,1:cnt3]
amax=max(strip)
}
if (cnt2 >= zdim[1,,]+1) {
count=zdim[1,,]+1
}
count=count+1
}
}
sortband_cube[x_cnt,y_cnt,]=sortband[1,1,]
#return(sortband_cube)
y_cnt=y_cnt+1
}
x_cnt=x_cnt+1
y_cnt=1
}
return(int(sortband_cube))
}
define sum_group_conc() {
#Sum mineral percentages from sma output, works with individual spectra and cubes
if ($ARGC == 0) {
printf ("\n")
printf (" $1 = smaout structure\n")
printf (" [$2]= .forcedlib from sma() output (must contain .group) \n")
printf (" Example: \n")
printf (" mineral_maps=sum_groups(final.endlib,final.conc,final.measured,final.nsamples,final.bb,final.forcedlib) \n")
printf (" \n")
printf (" Returns mineral group concentration maps with labels: \n")
printf (" .grouped_conc,.grouped_labels \n")
printf (" Also returns propagated group errors \n")
printf (" See also summary_sma_groups()--gives summary table of grouped abundances for a \n")
printf (" given x,y position. \n")
printf (" \n")
printf (" Groups are formed by scanning .group, and forming a new group index every \n")
printf (" time it encounters a unique group label. So, 'sheet-silicate' and 'sheet silicate' would \n")
printf (" be different groups. If you want to combine groups, like 'clays' and 'glass', you have \n")
printf (" to change both labels so that they match (such as 'clays + glass') and re-run sum_group_conc.\n")
return(0)
}
smaout=$1
end2_lib=smaout.endlib
lib_groups=end2_lib.group
min_conc=smaout.conc
endspc=translate(end2_lib.data,y,z)
meas_spc=smaout.measured
numsamp=smaout.nsamples
bbc=smaout.bb
low_chan=smaout.spectral_channel_low
high_chan=smaout.spectral_channel_high
# old method of input...kinda messy
# end2_lib=$1
# lib_groups=end2_lib.group
# min_conc=$2
# endspc=translate(end2_lib.data,y,z)
# meas_spc=$3
# numsamp=$4
# bbc=$5
# low_chan=$6
# high_chan=$7
if($ARGC==2) {
atlib=$2
lib_groups=cat(lib_groups,atlib.group,y)
atmspc=translate(atlib.data,y,z)
endspc=cat(endspc,atmspc,x)
}
meas_spc=meas_spc[,,low_chan:high_chan]
endspc=endspc[,low_chan:high_chan]
min_conc=cat(min_conc,bbc,z)
lib_groups=cat(lib_groups,"Black body",y)
endspc=cat(endspc,endspc[1]/endspc[1],x)
if(length(lib_groups)!=dim(min_conc)[3]) {
printf("Number of group labels must equal z-direction of .conc. Perhaps you have an atmospheric library to add? \n")
printf("Type sum_group_conc() for help \n")
return(0)
}
#first lets find the unique elements in the label array
glabels=unique(lib_groups)
numgroups=length(glabels)
group_index=create(dim(min_conc)[3],1,1,format=int,start=0,step=0)
for(j=1;j<=length(glabels);j+=1) {
group_index+=translate((lib_groups==glabels[,j])*j,y,x)
}
sumconc=create(dim(min_conc)[1],dim(min_conc)[2],numgroups,format=float,start=0.0,step=0)
grp_error=create(dim(min_conc)[1],dim(min_conc)[2],numgroups,format=float,start=0,step=0)
#Sum the groups and errors
for(k=1;k<=dim(min_conc)[1];k+=1) {
for(m=1;m<=dim(min_conc)[2];m+=1) {
newconc=create(1,1,1,format=float,start=0,step=0)
newlib=create(1,dim(endspc)[2],1,format=float,start=0,step=0)
newgrpindex=0
for(p=1;p<=dim(min_conc)[3];p+=1) {
if(min_conc[k,m,p]!=0.000000) {
newlib=cat(newlib,endspc[p],x)
newconc=cat(newconc,min_conc[k,m,p],y)
newgrpindex=cat(newgrpindex,group_index[p],x)
}
}
if(dim(newlib)[1]>1) {
newlib=newlib[2:]
newconc=newconc[,2:]
newgrpindex=newgrpindex[2:] # now have endlib and percentage vector for concentrations > 0.0
} else {
newlib=endspc
newconc=translate(min_conc[k,m],y,z)
}
# Contruct covariance matrix
em=translate(meas_spc[k,m,],x,z)
emT=translate(em,x,y)
newconcT=translate(newconc,x,y)
newlibT=translate(newlib,x,y)
term1=mxm(em,emT)
term2=mxm(newconcT,newlibT)
term3=mxm(term2,emT)
dimnewlib=dim(newlib)
term4=numsamp-dimnewlib[1]
term5=minvert(mxm(newlibT,newlib))
covm=((term1-term3)/term4)*(term5)
#Now, from those whittled down matrices, sum the concentrations and errors for each group
if(newgrpindex[1]!=0) {
for(i=1;i<=dim(newgrpindex)[1];i+=1) {
tmp_conc=create(1,1,1,format=float,start=0.0,step=0)
for(j=1;j<=dim(newconc)[2];j+=1) {
if(newgrpindex[j]==newgrpindex[i]) {
tmp_conc=float(cat(tmp_conc,float(newconc[,j]),y))
grp_error[k,m,(newgrpindex[i])]=sum(cat(grp_error[k,m,(newgrpindex[i])],float(covm[i,j]),x),x)
}
}
sumconc[k,m,(newgrpindex[i])]=sum(tmp_conc,y)
}
}
}
}
gp=struct()
gp.grouped_conc=sumconc
gp.grouped_error=sqrt(abs(grp_error))
gp.grouped_labels=glabels
return(gp)
}
define summary_sma(ylow,yhigh,xlow,xhigh,group,debug,output) {
#re-written by c.edwards 8/25/10
#combined summary_sma and summary_sma_group (automatically normalizes for bb, not atm)
if (($ARGC!=1 && $ARGC!=3) || ($ARGC==1 && HasValue(xlow)==0 && HasValue(ylow)==0)) {
printf (" Make summary table of results for one pixel from \n")
printf (" sma() output. \n")
printf (" $1=output structure from sma() \n")
printf (" $2=x-coordinate \n")
printf (" $3=y-coordinate \n")
printf (" xlow=low value for x-range\n")
printf (" xhigh=high value for x-range\n")
printf (" ylow=low value for y-range\n")
printf (" yhigh=high value for y-range\n")
printf (" NOTE: either a range or coordinate must be entered\n")
printf (" \n")
printf (" Note that only concentrations > 0.0001 are printed out \n\n")
printf (" Options:\n")
printf (" group=1 print the .grouped abundances if they are present\n")
printf (" output=1 return a structure with the sorted abundances, error, and endmember names\n")
printf (" \n")
return(null)
}
if(HasValue(debug)==0) { # Debug option added by S. Marshall, 2010-11-09
debug=0
}
inpu=$1
if(HasValue(group)==0) {
group=0
}
if($ARGC>1) {
xcor=$2
ycor=$3
if(xcor>dim(inpu.conc)[1]) xcor=dim(inpu.conc)[1]
if(ycor>dim(inpu.conc)[2]) ycor=dim(inpu.conc)[2]
}
#begin handling the range case by checking to make sure a valid range is input
if((HasValue(xlow) && HasValue(xhigh)==0) || (HasValue(xhigh) && HasValue(xlow)==0) && $ARGC==1) {
printf("You must enter both a xlow and xhigh for a complete range\n")
return(null)
}
if((HasValue(ylow) && HasValue(yhigh)==0) || (HasValue(yhigh) && HasValue(ylow)==0) && $ARGC==1) {
printf("You must enter both a ylow and yhigh for a complete range\n")
return(null)
}
#set defaults for both 3 arg case and 1 arg case
if(HasValue(ylow)==0 && HasValue(ycor)) {
ylow=ycor
} else if (HasValue(ylow)==0) {
ylow=1
}
if(HasValue(yhigh)==0 && HasValue(ycor)){
yhigh=ycor
} else if (HasValue(yhigh)==0) {
yhigh=1
}
if(HasValue(xlow)==0 && HasValue(xcor)) {
xlow=xcor
} else if(HasValue(xlow)==0) {
xlow=1
}
if(HasValue(xhigh)==0 && HasValue(xcor)) {
xhigh=xcor
} else if(HasValue(xhigh)==0) {
xhigh=1
}
#make sure we don't excede the bounds of the data cube
if(HasValue(ylow)) {
if(ylow<1) ylow=1
if(ylow>dim(inpu.conc)[2]) ylow=dim(inpu.conc)[2]
}
if(HasValue(yhigh)){
if(yhigh>dim(inpu.conc)[2]) yhigh=dim(inpu.conc)[2]
}
if(HasValue(xlow)){
if(xlow<1) xlow=1
if(xlow>dim(inpu.conc)[1]) xlow=dim(inpu.conc)[1]
}
if(HasValue(xhigh)){
if(xhigh>dim(inpu.conc)[1]) xhigh=dim(inpu.conc)[1]
}
labelflag=0
errorflag=0
#extract the correct elements for the structure so we have consitent variable names
if(group==1 && HasValue(get_struct(inpu,"grouped"))==1) {
#the grouped case
conc=inpu.grouped.grouped_conc[xlow:xhigh,ylow:yhigh]
if(HasValue(inpu.grouped.grouped_normconc)) {
norm=inpu.grouped.grouped_normconc[xlow:xhigh,ylow:yhigh]
} else {
norm=conc*0.
}
sort=inpu.grouped.grouped_sort_full[xlow:xhigh,ylow:yhigh]
label=inpu.grouped.grouped_labels
if(HasValue(get_struct(inpu.grouped,"grouped_bb"))) {
bb=inpu.grouped.grouped_bb[xlow:xhigh,ylow:yhigh]
} else {
bb=0.
}
if(HasValue(get_struct(inpu.grouped,"grouped_error"))) {
error=inpu.grouped.grouped_error[xlow:xhigh,ylow:yhigh]
bberror=inpu.grouped.grouped_bberror[xlow:xhigh,ylow:yhigh]
errorflag=1
}
} else {
#the full data case
conc=inpu.conc[xlow:xhigh,ylow:yhigh]
if(HasValue(inpu.normconc)) {
norm=inpu.normconc[xlow:xhigh,ylow:yhigh]
} else {
norm=conc*0.
}
sort=inpu.sort_full[xlow:xhigh,ylow:yhigh]
if(inpu.endlib.label[,1]!="none") {
label=inpu.endlib.label
}
if(length(inpu.forcedlib.label)!=1) {
if(inpu.forcedlib.label[,1]!="none") {
label=cat(label,inpu.forcedlib.label,y)
}
}
if(HasValue(get_struct(inpu,"bb"))) {
bb=inpu.bb[xlow:xhigh,ylow:yhigh]
} else {
bb=0.
}
if(HasValue(get_struct(inpu,"error"))) {
error=inpu.error[xlow:xhigh,ylow:yhigh]
bberror=inpu.bberror[xlow:xhigh,ylow:yhigh]
errorflag=1
}
}
rms=inpu.rms[xlow:xhigh,ylow:yhigh]
#if we don't have any labels
if(HasValue(label)==0) {
label=create(1,dim(conc)[3],1,start=1,step=1,format=int)
labelflag=1
}
for(i=1;i<=xhigh-xlow+1;i+=1) {
for(j=1;j<=yhigh-ylow+1;j+=1) {
xcor=i
ycor=j
cnt=1
if (debug == 1) {
printf("xlow = %d, ylow = %d, xcor = %d, ycor = %d\n", xlow, ylow, xcor, ycor)
}
printf ("\n")
printf ("*********************Summary for Pixel Coordinate (%i,%i) ******************* \n", xlow+xcor-1,ycor+ylow-1)
if(HasValue(get_struct(inpu,"label"))) {
if(ycor+ylow-1==1 && xcor+xhigh-1!=1) {
printf ("\t\t Label: %s\n\n",inpu.label[,xcor+xlow-1])
} else {
printf ("\t\t Label: %s\n\n",inpu.label[,ycor+ylow-1])
}
}
if(HasValue(get_struct(inpu,"sample_name"))) {
if(ycor+ylow-1==1 && xcor+xhigh-1!=1) {
printf ("\t\t Label: %s\n\n",inpu.sample_name[,xcor+xlow-1])
} else {
printf ("\t\t Label: %s\n\n",inpu.sample_name[,ycor+ylow-1])
}
}
if(group==1) {
printf (" Group Abundance(%%) Normalized for BB(%%) \n")
} else {
printf (" Endmember Abundance(%%) Normalized for BB(%%) \n")
}
printf("\n")
#setup the output structure
struct={}
if(labelflag==0) {
struct.label=text(1)
} else {
struct.label=0
}
if(errorflag==1) {
struct.error=0.
struct.normerror=0.
}
struct.conc=0.
struct.normconc=0.
struct.endlib=create(1,1,dim(inpu.endlib.data)[3],format=float,start=0,step=0)
if(dim(inpu.endlib.data)[2]!=1) {
inpu.endlib.data=translate(inpu.endlib.data,x,y)
}
while(sort[i,j,cnt]>0 && cnt<=dim(sort)[3]) {
idx=sort[i,j,cnt]
if(conc[,,idx] !=0) {
if(conc[,,idx] > 0.0001 || conc[,,idx] < -0.0001) {
if(labelflag==0 && errorflag==1) {
printf ("%27s %11.2f +/- %5.2f %13.2f +/- %5.2f\n", label[1:25,idx], conc[i,j,idx]*100, error[i,j,idx]*100, norm[i,j,idx], (error[i,j,idx]/sum(conc[i,j],axis=z))*100)
} else if (labelflag==1 && errorflag==1) {
printf ("%27i %11.2f +/- %5.2f %13.2f +/- %5.2f\n", label[,idx], conc[i,j,idx]*100, error[i,j,idx]*100, norm[i,j,idx], (error[i,j,idx]/sum(conc[i,j],axis=z))*100)
} else if (labelflag==0 && errorflag==0) {
printf ("%27s %11.2f %18.2f \n", label[1:25,idx], conc[i,j,idx]*100, norm[i,j,idx])
} else if (labelflag==1 && errorflag==0) {
printf ("%27i %11.2f %18.2f \n", label[,idx], conc[i,j,idx]*100, norm[i,j,idx])
}
#output the data in a useful format
if(labelflag==0) {
struct.label=cat(struct.label,label[1:25,idx],axis=y)
} else {
struct.label=cat(struct.label,label[,idx],axis=y)
}
struct.conc=cat(struct.conc,conc[i,j,idx]*100,axis=x)
struct.normconc=cat(struct.normconc,norm[i,j,idx],axis=x)
struct.endlib=cat(struct.endlib,inpu.endlib.data[idx],axis=x)
if(errorflag==1) {
struct.error=cat(struct.error,error[i,j,idx]*100,axis=x)
struct.normerror=cat(struct.normerror,float(error[i,j,idx]/sum(conc[i,j],axis=z)*100),axis=x)
}
}
}
cnt=cnt+1
}
printf ("___________________________________________________________________________\n")
printf (" Sum (BB included) = %5.2f\n", float(sum(conc[i,j],z)+bb[i,j])*100)
if(sum(norm[i,j],z)!=0) {
printf (" Sum (BB normalized) = %5.2f\n", float(sum(norm[i,j],z)))
}
struct.bbsum=float(sum(conc[i,j],z)+bb[i,j])*100
struct.bbabundance=float(bb[i,j]*100)
if(errorflag==1) {
printf (" Blackbody Abundance = %5.2f +/- %5.2f\n", float(bb[i,j]*100), float(bberror[i,j]*100))
struct.bberror=float(bberror[i,j]*100)
}
if(errorflag==0) {
printf (" Blackbody Abundance = %5.2f \n", float(bb[i,j]*100))
}
printf (" RMS Error(%%) = " + float(rms[i,j]*100) + " \n")
printf ("\n")
}
}
if(HasValue(output)) {
struct.conc=struct.conc[2:]
struct.normconc=struct.normconc[2:]
struct.label=struct.label[,2:]
struct.endlib=struct.endlib[2:]
if(errorflag==1) {
struct.error=struct.error[2:]
struct.normerror=struct.normerror[2:]
}
return(struct)
}
}
define summary_sma_group(ylow,yhigh,xlow,xhigh) {
if (($ARGC!=1 && $ARGC!=3) || ($ARGC==1 && HasValue(xlow)==0 && HasValue(ylow)==0)) {
printf ("\n Make summary table of abundances grouped by mineral class \n")
printf (" for one pixel from sma() output. \n")
printf (" $1=sma() output (must have .grouped) \n")
printf (" $2=x-coordinate \n")
printf (" $3=y-coordinate \n")
printf (" xlow=low value for x-range\n")
printf (" xhigh=high value for x-range\n")
printf (" ylow=low value for y-range\n")
printf (" yhigh=high value for y-range\n")
printf (" NOTE: either a range or coordinate must be entered\n")
printf (" \n")
return(null)
}
inpu=$1
if($ARGC>1) {
xcor=$2
ycor=$3
if(xcor>dim(inpu.conc)[1]) xcor=dim(inpu.conc)[1]
if(ycor>dim(inpu.conc)[2]) ycor=dim(inpu.conc)[2]
}
#begin handling the range case by checking to make sure a valid range is input
if((HasValue(xlow) && HasValue(xhigh)==0) || (HasValue(xhigh) && HasValue(xlow)==0) && $ARGC==1) {
printf("You must enter both a xlow and xhigh for a complete range\n")
return(null)
}
if((HasValue(ylow) && HasValue(yhigh)==0) || (HasValue(yhigh) && HasValue(ylow)==0) && $ARGC==1) {
printf("You must enter both a ylow and yhigh for a complete range\n")
return(null)
}
#set defaults for both 3 arg case and 1 arg case
if(HasValue(ylow)==0 && HasValue(ycor)) {
ylow=ycor
} else if (HasValue(ylow)==0) {
ylow=1
}
if(HasValue(yhigh)==0 && HasValue(ycor)){
yhigh=ycor
} else if (HasValue(yhigh)==0) {
yhigh=1
}
if(HasValue(xlow)==0 && HasValue(xcor)) {
xlow=xcor
} else if(HasValue(xlow)==0) {
xlow=1
}
if(HasValue(xhigh)==0 && HasValue(xcor)) {
xhigh=xcor
} else if(HasValue(xhigh)==0) {
xhigh=1
}
#make sure we don't excede the bounds of the data cube
if(HasValue(ylow)) {
if(ylow<1) ylow=1
if(ylow>dim(inpu.conc)[2]) ylow=dim(inpu.conc)[2]
}
if(HasValue(yhigh)){
if(yhigh>dim(inpu.conc)[2]) yhigh=dim(inpu.conc)[2]
}
if(HasValue(xlow)){
if(xlow<1) xlow=1
if(xlow>dim(inpu.conc)[1]) xlow=dim(inpu.conc)[1]
}
if(HasValue(xhigh)){
if(xhigh>dim(inpu.conc)[1]) xhigh=dim(inpu.conc)[1]
}
#we can now just call summary_sma with the group=1 option
summary_sma(inpu,xlow=xlow,xhigh=xhigh,ylow=ylow,yhigh=yhigh,group=1)
}
define tdb2struct(lab1, lab2) {
if ($ARGC == 0) {
printf("\nConvert tdb structure or file path to hdf structure\n")
printf(" usage: a = tdb2struct(tdb_file or tdb_file_path, [lab1=1][lab2=1])\n")
printf(" \n")
printf(" where:\n")
printf(" tdb_file = tdb format structure (load_tdb option distribute=1) of data, labels, and xaxes\nOR\n")
printf(" tdb_file_path = path to tdb format file of data, labels, and xaxes\n")
printf(" options:\n")
printf(" [lab=1] resample to xlab1 (923 pt lab spectra (221.7 to 1999.8) (default))\n")
printf(" [lab=2] resample to xlab2 (1971 pt lab spectra (200.5 to 3999.7))\n")
printf("\np. christensen 11/24/2010\n\n")
return(null)
}
struc = struct()
if(type($1)=="STRING") {
a=load_tdb($1,distribute=1)
} else if(type($1)=="STRUCT") {
a = $1
} else {
printf("$1 is either not a STRUCT or a PATH to a TDB file\n")
return(null)
}
n = length(a)
xlab1 = make_band(lab1=1)
xlab2 = make_band(lab2=1)
if(HasValue(lab2)) {
# resample to xlab2
data = resample(translate(a[1].data,x,z), translate(a[1].xaxis,x,z), xlab2).data
xaxis = xlab2
} else if(HasValue(lab1)) {
# resample to xlab1 - default
data = resample(translate(a[1].data,x,z), translate(a[1].xaxis,x,z), xlab1).data
xaxis = xlab1
} else {
# resample to xlab1 - default
data = resample(translate(a[1].data,x,z), translate(a[1].xaxis,x,z), xlab1).data
xaxis = xlab1
}
label = rm_quote(a[1].label)
for(i=2; i<=n; i+=1) {
if(HasValue(lab2)) {
# resample to xlab2
datan = resample(translate(a[i].data,x,z), translate(a[i].xaxis,x,z), xlab2).data
} else if(HasValue(lab1)) {
# resample to xlab1 - default
datan = resample(translate(a[i].data,x,z), translate(a[i].xaxis,x,z), xlab1).data
} else {
# resample to xlab1 - default
datan = resample(translate(a[i].data,x,z), translate(a[i].xaxis,x,z), xlab1).data
}
data = cat(data, datan, x)
label = cat(label, rm_quote(a[i].label), y)
}
struc.data = data
struc.label = label
struc.xaxis = xaxis
return(struc)
}
define fit_bb(type,tlow,thigh,steps,wave1,wave2,nn,instT,instdT,diffbb) {
if($ARGC<2) {
printf("\nGeneric Blackbody fitting routine\n")
printf("By default this uses a non-negative least squares fitting routine\n")
printf("to model a psudo-blackbody curve with many blackbodies\n")
printf("The output is based on the sma() output\n\n")
printf("\t$1 is the data that will be fit\n")
printf("\t$2 is the corresponding xaxis data in wavenumbers (1/cm)\n\n")
printf("Options:\n")
printf("\ttlow =\tthe low temperature range to fit (Default is 273.15)\n")
printf("\tthigh =\tthe high temperature range to fit (Default is 373.15)\n")
printf("\tsteps =\tthe number of steps between tlow and thigh (Default is 100)\n")
printf("\ttype =\t\"inst\" sets default values of tlow=223, thigh=323, steps=100\n")
printf("\ttype =\t\"cold\" sets default values of tlow=50, thigh=200, steps=100\n")
printf("\ttype =\t\"warm\" sets default values of tlow=150, thigh=350, steps=100\n")
printf("\ttype =\t\"hot\" sets default values of tlow=250, thigh=500, steps=100\n")
printf("\ttype =\t\"full\" sets default values of tlow=50, thigh=500, steps=100\n")
printf("\tinstT = sets the middle instrument temperature (default is (tlow+thigh)/2)\n")
printf("\tinstdT = sets the deviation around the instrument temperature (default is 20)\n")
printf("\twave1 =\tminimum wavenumber for modeling (Default is min(xaxis))\n")
printf("\twave2 =\tmaximum wavenumber for modeling (Default is max(xaxis))\n")
printf("\tnn =\tuse the non-negative fitting routine (Default is 1)\n")
printf("\tdiffbb = models the differencd of two blackbodies using the same temperature ranges set by \"type\"\n")
printf("\t\tIt creates a range of blackbody temperature differences from thigh to tlow\n\n")
printf("c.edwards 5-5-2011\n\n")
return(null)
}
data=$1
xaxis=$2
maxsteps=dim(xaxis)[3]-1
#set the temperature ranges based on the desired fits
if(HasValue(type)) {
if(type=="inst") {
if(HasValue(tlow)==0) tlow=223
if(HasValue(thigh)==0) thigh=323
} else if(type=="cold") {
if(HasValue(tlow)==0) tlow=50
if(HasValue(thigh)==0) thigh=200
} else if(type=="warm") {
if(HasValue(tlow)==0) tlow=150
if(HasValue(thigh)==0) thigh=350
} else if(type=="hot") {
if(HasValue(tlow)==0) tlow=250
if(HasValue(thigh)==0) thigh=500
} else if(type=="full") {
if(HasValue(tlow)==0) tlow=50
if(HasValue(thigh)==0) thigh=500
} else {
printf("No valid location for defaults was entered\n")
printf("Please try again\n")
return(null)
}
}
#if we have a standard dT to be used, set it instead of a range
if(HasValue(instT)) {
if(HasValue(instdT)==0) {
instdT=20
}
if(HasValue(tlow)==0) tlow=instT-instdT
if(HasValue(thigh)==0) thigh=instT+instdT
}
#use nnls by default as this yeilds superior results in this case
#and in many cases, it is the only one that works
if(HasValue(nn)==0) {
nn=1
}
#set the wavelength ranges for fitting
if(HasValue(wave1)==0 && HasValue(wave2)==0) {
wave1=min(xaxis)
wave2=max(xaxis)
}
#finally if we don't have any other options, set the range
if(HasValue(tlow)==0) {
tlow=273.15
}
if(HasValue(thigh)==0) {
thigh=373.15
}
#set the number of steps
if(HasValue(steps)==0) {
steps=100
}
#check to see if we have too many steps
if(steps>maxsteps) {
printf("Error: more steps than allowed\n")
printf("Resetting to Xaxis length-1\n")
steps=maxsteps-1
}
#now we see if we are modeling bb differences or straight bbs
if(HasValue(diffbb)==0) {
#setup the bb range for modeling
bb={}
bb.data=create(steps,1,dim(xaxis)[3],start=0,step=0,format=float,org=bsq)
bb.temp=create(steps,1,1,start=0,step=0,format=float,org=bsq)
bb.label=text(steps)
bb.xaxis=xaxis
#temperature ranges based on the number of steps and tlow/thigh
for(i=1;i<=steps;i+=1) {
temp=tlow+(i-1)*((thigh-tlow)/float(steps))
bb.data[i]=bbr(xaxis,temp)
bb.label[,i]=sprintf("%.2fK Blackbody",temp)
bb.temp[i]=temp
}
} else {
#this is a little more complicated
#we take the sqrt of the number of steps to do an xy array difference of the
#total number of steps
stepssq=int(ceil(sqrt(steps)))
#check to see if we have too many steps
if(stepssq^2>maxsteps) {
stepssq=int(floor(sqrt(maxsteps-1)))
}
steps=stepssq^2
#same setup as before
bb={}
bb.data=create(steps,1,dim(xaxis)[3],start=0,step=0,format=float,org=bsq)
bb.temp=create(steps,2,1,start=0,step=0,format=float,org=bsq)
bb.xaxis=xaxis
#make the 2-D array of temperature differences
count=1
for(j=1;j<=stepssq;j+=1) {
for(i=1;i<=stepssq;i+=1) {
temp1=tlow+(i-1)*((thigh-tlow)/float(stepssq))
temp2=thigh-(j-1)*((thigh-tlow)/float(stepssq))
bb.data[count]=bbr(xaxis,temp2)-bbr(xaxis,temp1)
bb.temp[count]=cat(temp2,temp1,axis=y)
count++
}
}
#extract the values that have zero difference (this makes sma unhappy)
mask=sum(bb.data,axis=z)!=0
bb.data=extract(bb.data,mask)
bb.temp=extract(bb.temp,mask)
newsteps=dim(bb.temp)[1]
bb.label=text(newsteps)
#fill out the label after we remove the zero differences
for(i=1;i<=dim(bb.temp)[1];i+=1) {
bb.label[,i]=sprintf("%.2fK Blackbody - %.2fK Blackbody",bb.temp[i,1],bb.temp[i,2])
}
}
#model the data
out=sma(data,bb,bb=0,group=0,wave1=wave1,wave2=wave2,nn=nn,noerr=1)
#setup some output stats
if(HasValue(diffbb)==0) {
#setup the output
numpos=int(sum(out.sort!=0))
shorttemp=shortconc=create(numpos,1,1,start=0,step=0,format=float,org=bsq)
shortrad=create(numpos,1,dim(xaxis)[3],start=0,step=0,format=float,org=bsq)
#get only data that are non-zero
for(i=1;i<=numpos;i+=1) {
shorttemp[i]=bb.temp[out.sort[,,i]]
shortconc[i]=out.conc[,,out.sort[,,i]]
shortrad[i]=bb.data[out.sort[,,i]]
}
#here we calculate the normalized concentrations based on power and temperature
global(con)
shortpow=con.sigma*(shorttemp^4)
scaleconc=shortconc/shortpow
normconc=scaleconc/sum(scaleconc)
#output the data including a weighted average temperature
stats={}
stats.pow=shortpow
stats.normconc=normconc
stats.conc=shortconc
stats.temp=shorttemp
stats.rad=shortrad
stats.avgtemp=sum(stats.normconc*stats.temp)
} else {
#setup the additional output for diffbb
numpos=int(sum(out.sort!=0))
shorttemp=create(numpos,2,1,start=0,step=0,format=float,org=bsq)
shortconc=create(numpos,1,1,start=0,step=0,format=float,org=bsq)
#get the non-zero data
for(i=1;i<=numpos;i+=1) {
shorttemp[i]=bb.temp[out.sort[,,i]]
shortconc[i]=out.conc[,,out.sort[,,i]]
}
#calcualate a normalized concentration
normconc=shortconc/sum(shortconc)
#output the data into the stats structure
stats={}
stats.temp=shorttemp
stats.conc=shortconc
stats.normconc=normconc
}
out.stats=stats
return(out)
}
define sa(deg) {
if($ARGC == 0) {
printf("solid angle for phi from 0 to phi and theta from theta1 to theta2\n")
printf("usage: a = sa(phi, theta1(rad), theta2(rad) [deg])\n")
printf("Default is in radians; use deg=1 to input in degrees\n")
return(0)
}
# p. christensen 3/10
if(HasValue(deg)) {
phi = d2r($1)
t1 = d2r($2)
t2 = d2r($3)
printf("phi=%.2f t1= %.2f t2=%.2f\n", phi, t1, t2)
} else {
phi = $1
t1 = $2
t2 = $3
}
omega = phi * ( cos(t1) - cos(t2)); # steradians
return(omega)
}