a strange one

[Pixie Songbook]

This is the program I have, it takes data from a file with two seperate columns of data and attempts to find the best model for the data using a chi-squared minimisation procedure. The odd thing is if I have the program in the form below, it quits after only 6 iterations saying chi-squared has been minimised. If however I edit out line 119 and include line 123 instead the program continues doing as manby iterations as I tell it too. Can anyone tell me;
a) why the program stops so soon for the first scenario?
b) why the difference between the scenarios? :confused:
thanks
pixie
--------------------------------------------------------------------------------------------------------------------
#include <iostream.h>
#include <fstream.h>
#include <math.h>

int main()
{
int m, j;
float lambda[9096], L, I[9096]; // note: lambda[9096] is a lie, 4096 not used as it result in calculation error
ifstream read("C:\\Documents and Settings\\Administrator\\My Documents\\My Honours\\qso.doc");
ofstream prt("C:\\Documents and Settings\\Administrator\\My documents\\My Honours\\average wavelength, standard deviation and maximum intenesity value for qso data.doc", ios::trunc);
ofstream prt2("C:\\Documents and Settings\\Administrator\\My documents\\My Honours\\intensity.doc", ios::trunc);
ofstream prt3("C:\\Documents and Settings\\Administrator\\My documents\\My Honours\\wavelength.doc", ios::trunc);
if(!read)
{
cout << "input file failed to open";
system ("pause");
exit(1);
}
if(!prt)
{
cout << "output file 1 failed to open";
system ("pause");
exit(1);
}
if(!prt2)
{
cout << "output file 2 failed to open";
system ("pause");
exit(1);
}
if(!prt3)
{
cout << "output file 3 failed to open";
system ("pause");
exit(1);
}

cout << "the value of m is "; //use 2048 to calculate average over all wavelengths
cin >> m;

for (j = 0; j < m; j++) //extracting values from file
{
read >> L;
if(L - 2000 > 0) //this is to ensure only wavelengths are equated to lambda[j]
lambda[j] = L;
else
I[j] = L;

read >> L;
if(!(L - 2000 > 0)) //this is to ensure only intensities are equated to I[j]
I[j] = L;
else
lambda[j] = L;

cout << "lambda[" << j + 1 << "] = " << lambda[j] << " "; //(j + 1) so zero lambda[0] isn't first term
//prt << "lambda[" << j + 1 << "] = " << lambda[j] << endl;
prt3 << lambda[j] << endl;
cout << "I[" << j + 1 << "] = " << I[j] << endl; //(j + 1) so zero I[0] isn't first term, I[1] is
//prt << "I[" << j + 1 << "] = " << I[j] << endl;
prt2 << I[j] << endl;
}
read.close();
prt2.close();
prt3.close();

float average; //calculating average
average = 0;
for (j = 0; j < m; j++)
average = average + (lambda[j]/m);
cout << endl << "average wavelength = " << average << endl;
prt << "average wavelength = " << average << endl;

double Imax; //determining maximum intenisity value
Imax = 0;
for (j = 0; j < m; j++)
{
if(I[j] - Imax > 0)
Imax = I[j];
else if(I[j] - Imax <= 0)
Imax = Imax;
}
cout << "maximum intensity value = " << Imax << endl;
prt << "maximum intensity value = " << Imax << endl;

float sigma, S; //calculating standard deviation
sigma = 0;
for (j = 0; j < m; j++)
sigma = sigma + (((average - lambda[j])*(average - lambda[j]))/m);
S = sqrt(sigma);
cout << "standard deviation = " << S << endl;
prt << "standard deviation = " << S << endl;

ofstream prt4("C:\\Documents and Settings\\Administrator\\My documents\\My Honours\\gaussian.doc", ios::trunc);
if(!prt4)
{
cout << "output file 4 failed to open";
system ("pause");
exit(1);
}
ofstream intensity("C:\\Documents and Settings\\Administrator\\My documents\\My Honours\\intensity(calculated).doc", ios::trunc);
if(!intensity)
{
cout << "output file 'intensity' failed to open";
system ("pause");
exit(1);
}

double K_old, K_new, K_min; //estimates of K
int n_max, n; //where n_max = number of iterations & n = running index

cout << endl << "enter your initial guess for K ... " << endl;
cin >> K_old; //guess about 450?

cout << endl << "enter the maximum number of iterations ... " << endl;
cin >> n_max;

double chi_min;
chi_min = 2000;
double chi[n_max], chi_prime[n_max];
for(n = 0; n <= n_max ; n++)
{
double Chi, Chi_prime;
//double chi[n], chi_prime[n];
double X, i[9096];
Chi = 0;
Chi_prime = 0;
X = 2 * S * S * -1; //needed for gaussian form below
for (j = 0; j < m; j++)
{
i[j] = Imax * exp((lambda[j] - average)/X) - K_old ; //use K = ? here
Chi = Chi + pow(((I[j] - i[j])/S), 2);
Chi_prime = Chi_prime + 2*((I[j] - i[j])/S)*(1/S); //derivative of chi with respect to K
}
chi[n] = Chi;
chi_prime[n] = Chi_prime;
cout << "chi squared = " << chi[n] << " ";
prt4 << "chi squared = " << chi[n] << " ";
prt4 << "derivative = " << chi_prime[n] << " ";
K_new = K_old - (chi[n]/chi_prime[n]);
cout << "K = " << K_new << endl;
if (0 < K_new && K_new < 1368)
{
if (chi_min - chi [n] > 0)
{
chi_min = chi[n];
K_min = K_new;
prt4 << "K = " << K_new << endl << endl;
}
else {
chi_min = chi_min;
prt4 << "chi not minimal" << endl << endl;
}
}
else {
chi_min = chi_min;
prt4 << "K < 0 or K > 1368" << endl << endl;
}

if (n > 2)
{
if ((chi[n] = chi_min) && (fabs(chi[n - 1] - chi_min) < 0.1) && (fabs(chi[n - 2] - chi_min) < 0.1))
{
cout << "minimised chi = " << chi_min << " " << "corresponding K value = " << K_min << endl << endl;
prt << "minimised chi = " << chi_min << " " << "corresponding K value = " << K_min << endl;
cout << "chi squared minimised before " << n_max << " iterations completed" << endl;
prt4 << "chi squared minimised before " << n_max << " iterations completed" << endl;
for (j = 0; j < m; j++)
{
//cout << "i[" << j + 1 << "] = " << i[j] << endl;
//intensity << "i[" << j + 1 << "] = " << i[j] << endl;
intensity << i[j] << endl;
}
system ("pause");
exit(1);
}
}
K_old = K_new;
}
cout << "minimised chi = " << chi_min << " " << "corresponding K value = " << K_min << endl;
prt << "minimised chi = " << chi_min << " " << "corresponding K value = " << K_min << endl;
prt.close();
prt4.close();
system ("pause");
}