Error Characterization of the Alpha Residuals Emi ssivity Extraction Technique

Michael C. Baglivio

/****************************************************************************

TITLE: alpha_resid.c

BY: Mike Baglivio

CLASS: Senior Research

PURPOSE: To characterize error associated with Weins Approximation

in the alpha residuals T/e separation technique.

*****************************************************************************/

#include<iostream.h>

#include<fstream.h>

#include<stdlib.h>

#include<math.h>

#define Y_SIZE 5

#define X_SIZE 5

#define NUM_BANDS 128

const double Pi = 3.141592654;

const double h = 6.634E-34;

const double k = 1.38054E-23;

const double c = 3.0E+8;

const double con1 = 2*Pi*h*c*c;

const double con2 = c*h/k;

const int n = 128;

const long int input_data_points = 35350;

int main(int argc, char *argv[])

{

//PROMPTS USER FOR NECCESSARY INPUT

int sub_sample;

char outfile[20];

int width;

int temp = 300;

cout<<"Use ODD numbers for sample points and response width."<<endl;

cout<<"How many sample points per channel?> ";

cin>>sub_sample;

cout<<"How wide would you like the spectral response curves to be?> ";

cin>>width;

cout<<"Output file name> ";

cin>>outfile;

//DEFINES START WAVELENGTH OF SAMPLE_WAVE ARRAY

long int i;

int j;

int q=0;

double wave_incr = .04e-6;

double init_w_val = 7.80e-6;

double w_val = 7.80e-6;

//DEFINES SUBSAMPLING AND INTERVAL FOR SAMPLE POINTS IN EACH CHANNEL

int sample_points = sub_sample*n + (width-sub_sample);

double sub_sample_incr = wave_incr/sub_sample;

double *sample_wave = NULL;

sample_wave = new double[ sample_points ];

//CREATES OUTPUT FILE

ofstream fout(outfile);

//READS INPUT IMAGE FILE

int x, y;

static float spectrum[ NUM_BANDS ];

/* open the image file */

ifstream input_file( "test.img" );

if(!input_file){

cerr<<"Image file cannot be opened" <<endl;

exit(1);

}

for( y = 0; y <1/*Y_SIZE*/; y++ ) {

for( x = 0; x <1 /*X_SIZE*/; x++ ) {

input_file.read( spectrum, sizeof( float ) * NUM_BANDS );

//input_file.close();

//LOOP FOR LIBRARY CALCULATION

int file=0;

for(file=1; file<argc; file++){

//fout<<argv[file]<<endl;

q=0;

double wave_incr = .04e-6;

double init_w_val = 7.80e-6;

double w_val = 7.80e-6;

//PRODUCES GAUSSIAN CURVE TO MODEL RESPONSIVITY

double *gauss = NULL;

gauss = new double[ width ];

double g=0;

for(i=0; i<width; i++){

g=pow(Pi*(i-(width-1)/2),2);

gauss[i]=exp(-g/100);

}

//CALCULATION LOOP

int *center = NULL;

center = new int[ n ];

int a;

for(i=0; i<n; i++){

center[i]=(i*sub_sample)+(sub_sample-1)/2; //center of band

for(j=0; j<width; j++){

a=center[i]-(width-1)/2+j;

if(a<0)

a=0;

if(a>sample_points)

a=sample_points;

}

}

//MAKES CALCULATIONS TO COMPUTE ALPHA RESIDUALS OF INPUT DATA

double *X_one = NULL;

X_one = new double[ n ];

double bandwidth[n]; //wave[n] is the bandwidth array

init_w_val=7.80e-6;

double sum_X_one=0.0;

for(i=0; i<n; i++){

bandwidth[i]=init_w_val;

init_w_val = init_w_val + wave_incr;

X_one[i]=bandwidth[i]*log(spectrum[i]);

sum_X_one+=X_one[i];

}

double X_one_bar=sum_X_one/n;

double *alpha_one = NULL;

alpha_one = new double[ n ];

for(i=0; i<n; i++){

alpha_one[i]=X_one[i]-X_one_bar;

}

//THIS SECTION COMPUTES ALPHA RESIDUALS OF THE LIBRARY SPECTRA

//READS IN LIBRARY EMISSIVITY FILE

double *library = NULL;

library = new double[ input_data_points ];

double *lib_wave = NULL;

lib_wave = new double[ input_data_points ];

double *actual_library = NULL;

actual_library = new double[sample_points];

double *actual_lib_wave = NULL;

actual_lib_wave = new double[sample_points];

ifstream fin(argv[file]);

if(!fin){

cerr<<"File cannot be opened" <<endl;

exit(1);

}

i=0;

while(fin >>lib_wave[i]>>library[i]){

i++;

}

fin.close();

//CONVERTS LIB_WAVE[I] INTO MICRONS

for(i=0; i<input_data_points; i++){

lib_wave[i] = lib_wave[i]*10e-6;

}

//EXTRACTS OUT NECESSARY WAVELENGTH AND EMISSIVITY VALUES FROM INPUT DATA

for(i=0; i<sample_points; i++){

sample_wave[i]=w_val;

w_val = w_val + sub_sample_incr;

for(j=0; j<input_data_points; j++){

if(lib_wave[j]>=sample_wave[i] && lib_wave[j]<sample_wave[i]+0.0001){

actual_lib_wave[i] = lib_wave[j];

actual_library[i] = library[j];

}

}

}

//COMPUTES EFFECTIVE EMISSIVITY

double e_num;

double e_denom;

double *eff_emmiss = NULL;

eff_emmiss = new double[ n ];

double *X_two = NULL;

X_two = new double[ n ];

double sum_X_two=0.0;

int z=0;

for(i=0; i<n; i++){

e_num=e_denom=0.0;

for(j=0; j<width; j++){

z=center[i]-(width-1)/2+j;

if(z<0)

z=0;

if(z>sample_points)

z=sample_points;

e_num+=actual_library[z]*gauss[j]*wave_incr;

e_denom+=gauss[j]*wave_incr;

}

eff_emmiss[i]=e_num/e_denom;

X_two[i]=bandwidth[i]*log(eff_emmiss[i])+bandwidth[i]*log(con1)-5*bandwidth[i]*log(bandwidth[i])-bandwidth[i]*log(Pi)-con2/temp;

sum_X_two+=X_two[i];

//cout<<eff_emmiss[i]<<endl;

}

double X_two_bar=sum_X_two/n;

double *alpha_two = NULL;

alpha_two = new double[ n ];

for(i=0; i<n; i++){

alpha_two[i]=X_two[i]-X_two_bar;

//cout<<alpha_two[i]<<endl;

}

//COMPUTES DIFFERENCE OF INPUT AND LIBRARY SPECTRA

double difference;

double sum_diff=0.0;

double spec_diff=0.0;

for(i=0; i<n; i++){

sum_diff+=(alpha_two[i]-alpha_one[i])/alpha_two[i];

spec_diff=alpha_two[i]-alpha_one[i];

fout<<spec_diff<<endl;

}

difference=sum_diff/n;

} //END OF LIBRARY LOOP

//cout<<min<<endl;

}

}

return(0);

}