// Particle swarm optimization for measured aerosol optical depths at 13 wavelengths:
// 412, 440, 463, 479, 500, 520, 556, 610, 675, 750, 778, 870, and 1020 nm
// (If using 6 wavelengths consider only 440, 555, 674, 778, 869, 1019 with new solar spectrophotometer)
// version 4 adds scouts
// added bimodal lognormal distribution 15 Dec 2015
// add option to fix alpha or beta parameters 17 Dec 2015
// update 12/29/2015 to include testing for parameter fixed and only updating those that vary - affects init_swarm, init_particle, and update_swarm
// 12/29/2015 changed global parameter stop condition from 1.0E-10 to 1.0E-5
// c++ includes
#include <iostream>
#include <fstream>
#include <cstring>
#include <string>
#include <sstream>
#include <new>
#include <cmath>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <cstdio>

using namespace std;	// required for c++ file routines

// Put #defines here
#define channels 13 // number of wavelength channels to use for optimization
#define particles 200 // number of particles in the swarm
#define PI 3.1415192654		// the constant PI
#define CONV_FACT 1.0E-4	/* convert microns to cm */
#define CONV_FACT1 1.0E-3	// convert nanometers to micrometers
#define CONV_FACT2 1.0E-7	/* convert nanometers to cm */
#define c1 1.8 // acceleration coefficient 1
#define c2 1.8 // acceleration coefficient 2
#define pop_max 10 // define the maximum number of populations to get parameters for
#define gen_max 5000 // define the maximum number of generations to go through
#define N_tmax 100 // define the stop condition for number of trials with no significant change
#define scouts 10 // define the number of scouts used for additional global minimum checking

// required definition of matrix at top of source code
typedef vector<vector<double> > matrix; // define matrix object

// function prototype definitions to be placed before main() source code
// generates a random value between low and high numbers using log10 space
double random(double low, double high);
// calculates aod for each particle in the swarm -- call with calc_aod(swarm[x], mie[y], aodc[x][y])
void calc_aod(vector<double> &part, vector<double> &Qext, double &tauc, int dist);
// initialize mie extinction matrix by filling it with values from files
void initQ(vector<double> &Qext, string fname);
void init_mie(matrix &mie);
// calculate the reduced chi^2 value for each particle in the swarm
void calc_red_chi2(double aodm[], matrix &aodc, double f[], int dist);
// initialize the swarm by passing by passing range limits for parameters
void init_swarm(matrix &swarm, double param_min[], double param_max[], vector<double> &param, int num, int param_fix);
// generates a random particle using log10 space for parameter ranges - called by init_swarm()
void init_particle(vector<double> &part, double param_min[], double param_max[], vector<double> &param, int param_fix);
// find global reduced chi^2 minimum
int find_global_min(double f[], double Pg[], matrix &swarm, int num);
// find local reduced chi^2 minimum for each particle
void find_local_min(double f[], double f0[], matrix &Pl, matrix &swarm);
void update_swarm(matrix &swarm, matrix &Pl, double Pg[], double param_min[], double param_max[], vector<double> &param, int param_fix);


int main()
{
	bool d_global_small=false;	// boolean variable to test for significant changes in the global parameters
	int i,j; // generic indices for loops
	int population=1; // index for number of populations
	int dist; // type of distribution to fit to
	int i_g,gen=0; // indices for global minimum and generation number
	int i_j; // index for best extra particle's reduced chi^2
	int N_t=0; // numbers of trials with no significant change in parameters used for stop condition
	int param_fix=0; // holds info as to whether any parameters are to be fixed
	float day;
	string mystr; // holds typed in filenames
	char *fname1,*fname2,*fname3,dataline[500]; // holds input and output filenames and datalines from input file
	double aodm[channels]; // holds measured aerosol optical depths (AODs) for a particular date & time
	double f0[particles]; // initial reduced chi^2 vector
	double f[particles]; // current reduced chi^2 vector
	matrix swarm(particles,vector<double>(6,0.0)); // matrix of particles each with 6 dimensions N0, a, b, N1, a1, b1
	matrix aodc(particles,vector<double>(channels,0.0)); // matrix of calculated AODs for each particle at each wavelength
	matrix Pl(particles,vector<double>(6,0.0)); // matrix of local minimum adjustments
	double Pg[6]; // vector for global minimum adjustments for all particles
	double Pg0[6] = {1.0E12,0.1,0.307,1.0E11,1.0,0.307}; // initial global minima for checking stop conditions
	matrix Pj(scouts,vector<double>(6,0.0)); // vector for randomly generated particle so we don't get stuck in a local minimum - use matrix object to be compatible with swarm for calls to init_swarm
	double Pj0[6]; // holds the global minima for scouts
//	vector<double> aodj(channels,0.0);
	matrix aodj(scouts, vector<double>(channels,0.0)); //vector for extra particle's aods for global minimum searching improvement
	double fj[scouts]; // reduced chi^2 for extra randomly generated particles
//	float r[9961]; // vector of radii from 0.040-10.000 um in 0.001 increments
	matrix mie(channels,vector<double>(9961,0.0)); // matrix of extinction coefficients for each wavelength at each radius
	vector<double> param(6,0.0); // holds initial set of N0, alpha, and beta in order to get initial reduced chi^2 functions
//	double Nmin, Nmax, amin, amax, bmin, bmax; // range for parameters: N0, alpha, and beta
//	replace above line with arrays
	double param_min[6], param_max[6]; // minima and maxima for N0, alpha, and beta
	double red_chi2_goal; // relative eroor between calculated and measured aods
	fstream aodfile,outfile,outfile1;
	
	srand(time(NULL)); // seed the random number algorithm
	init_mie(mie);

	cout << "Enter measured AOD file to process (e.g., bimodal1.txt): ";
	getline(cin,mystr);
	fname1 = new char [mystr.size()+1];	/* c_string method */
	strcpy(fname1,mystr.c_str());

	cout << endl << "Enter Distribution type to fit:" << endl << "1 = Junge" << endl << "2 = Gamma" << endl << "3 = Lognormal" << endl << "4 = Bimodal Lognormal" << endl;
	cin >> dist;
	if(dist==1){red_chi2_goal=1.0E-6; cout << "Junge selected" << endl;} // Junge distribution
	else if(dist==2){red_chi2_goal=1.0E-5; cout << "Gamma selected" << endl;} // Gamma distribution
	else if(dist==3){red_chi2_goal=1.0E-5; cout << "Lognormal selected" << endl;} // Lognormal distribution
	else if(dist==4){red_chi2_goal=1.0E-3; cout << "Bimodal Lognormal selected" << endl;} // Bimodal Lognormal distribution
	else{cout << "Bad distribution!" << endl; return 0;}

	cout << "Enter output file name for best number distribution parameters N0, alpha, beta, N1, alpha1, beta1 (e.g., results2bimodal.txt): ";
	cin >> mystr;
	fname2 = new char [mystr.size()+1];	/* c_string method */
	strcpy(fname2,mystr.c_str());
	
	cout << "Enter output file name summary results: ";
	cin >> mystr;
	fname3 = new char [mystr.size()+1];	/* c_string method */
	strcpy(fname3,mystr.c_str());
	
	aodfile.open(fname1);	// open input file of aod values 
	if (!aodfile.is_open()){
		cout << "\nUnable to open file " << fname1 << "\n";
		return 0;
	}
	outfile.open(fname2,ios::app | ios::out);	// open output file for number distribution parameters N0, alpha, beta 
	if (!outfile.is_open()){
		cout << "\nUnable to open file " << fname2 << "\n";
		return 0;
	}
	outfile << "N0, alpha, beta, N1, alpha1, beta1" << endl;		// write header line to output file

	outfile1.open(fname3,ios::app | ios::out);	// open output file for summary results
	if (!outfile1.is_open()){
		cout << "\nUnable to open file " << fname3 << "\n";
		return 0;
	}
	outfile1 << "red. chi^2, N0, alpha, beta, N1, alpha1, beta1" << endl;		// write header line to output file

/*	This kept giving segmentation faults, so replaced with stringstream 
	while(!aodfile.eof()){
		getline(aodfile,mystr);			// get a line of data
		strcpy(dataline,mystr.c_str());		// extract data fields
		sscanf(dataline,"%f %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",&day,&aodm[0],&aodm[1],&aodm[2],&aodm[3],&aodm[4],&aodm[5],&aodm[6],&aodm[7],&aodm[8],&aodm[9],aodm[10],&aodm[11],&aodm[12]);		
	}
MAKE SURE aod files don't have commas, only spaces separating the values
*/
	while(getline(aodfile,mystr)){
		cout << mystr << endl;
		stringstream(mystr) >> day >> aodm[0] >> aodm[1] >> aodm[2] >> aodm[3] >> aodm[4] >> aodm[5] >> aodm[6] >> aodm[7] >> aodm[8] >> aodm[9] >> aodm[10] >> aodm[11] >> aodm[12];
	}

	for(i=0;i<channels;i++){
		cout << "aodm["<<i<<"] = " << aodm[i] << endl;;
	}

	cout << "closing file " << fname1 << endl;
	aodfile.close();
// ask for parameter ranges later, for now start with defaults
/*	printf("Input the range for N0 (e.g., 1.0E6,1.0E12): ");
	scanf("%lf,%lf", &Nmin,&Nmax);
	printf("Input the range for alpha (e.g., 0.01,1.0): ");
	scanf("%lf,%lf", &amin,&amax);
	printf("Input the range for beta (e.g., 0.1,1.0): ");
	scanf("%lf,%lf", &bmin,&bmax);
*/
	switch(dist){
		case 1:	// Junge distribution
		param[0]=1.0E6; // initial value for N0
		param_min[0]=1.0E3;
		param_max[0]=1.0E9;
		param[1]=3.0; // initial value for alpha (a)
		param_min[1]=2.0;
		param_max[1]=5.0;
		param[2]=0.30; // initial value for beta (b)
		param_min[2]=0.1;
		param_max[2]=0.5;
		param[3]=1.0E6; // initial value for N1
		param_min[3]=1.0E4;
		param_max[3]=1.0E8;
		param[4]=3.0; // initial value for alpha1 (a1)
		param_min[4]=2.0;
		param_max[4]=5.0;
		param[5]=0.30; // initial value for beta1 (b1)
		param_min[5]=0.1;
		param_max[5]=0.5;
		break;
		case 2:	// modified gamma distribution
		param[0]=1.0E9; // initial value for N0
		param_min[0]=1.0E6;
		param_max[0]=1.0E12;
		param[1]=0.1; // initial value for alpha (a)
		param_min[1]=0.05;
		param_max[1]=0.50;
		param[2]=0.35; // initial value for beta (b)
		param_min[2]=0.1;
		param_max[2]=0.5;
		param[3]=1.0E9; // initial value for N1
		param_min[3]=1.0E6;
		param_max[3]=1.0E12;
		param[4]=0.1; // initial value for alpha1 (a1)
		param_min[4]=0.05;
		param_max[4]=0.50;
		param[5]=0.35; // initial value for beta1 (b1)
		param_min[5]=0.1;
		param_max[5]=0.5;
		break;
		case 3:	//lognormal distribution
		param[0]=1.0E6; // initial value for N0
		param_min[0]=1.0E4;
		param_max[0]=1.0E10;
		param[1]=0.5; // initial value for alpha (a)
		param_min[1]=0.05;
		param_max[1]=1.0;
		param[2]=0.307; // initial value for beta (b)
		param_min[2]=0.1;
		param_max[2]=0.5;
		param[3]=1.0E6; // initial value for N1
		param_min[3]=1.0E4;
		param_max[3]=1.0E10;
		param[4]=0.5; // initial value for alpha1 (a1)
		param_min[4]=0.05;
		param_max[4]=1.0;
		param[5]=0.307; // initial value for beta1 (b1)
		param_min[5]=0.1;
		param_max[5]=0.5;
		break;
		default:	// default to bimodal lognormal distribution
		param[0]=1.0E7; // initial value for N0
		param_min[0]=1.0E5;
		param_max[0]=1.0E11;
		param[1]=0.1; // initial value for alpha (a)
		param_min[1]=0.05;
		param_max[1]=0.50;
		param[2]=0.307; // initial value for beta (b)
		param_min[2]=0.1;
		param_max[2]=0.5;
		param[3]=1.0E5; // initial value for N1
		param_min[3]=1.0E3;
		param_max[3]=1.0E9;
		param[4]=2.0; // initial value for alpha1 (a1)
		param_min[4]=0.5;
		param_max[4]=5.0;
		param[5]=0.307; // initial value for beta1 (b1)
		param_min[5]=0.1;
		param_max[5]=0.5;
	}

// Prompt user if parameters alpha and/or beta should be held constant
	cout << "Hold alpha or beta constant?" << endl;
	cout << "0: allow both alpha and beta to vary" << endl;
	cout << "1: hold alpha constant" << endl;
	cout << "2: hold beta constant" << endl;
	cout << "3: hold both alpha and beta constant" << endl;
	cout << "4: BIMODAL ONLY - hold alpha0 and beta0 constant; let alpha1 and beta1 vary" << endl;
	cin >> param_fix;
	switch (param_fix)
	{
		case 0:
			cout << "both alpha and beta will vary" << endl;
			outfile << "both alpha and beta will vary" << endl;
			break;
		case 1:
			param_min[1]=param[1];
			param_max[1]=param[1];
			param_min[4]=param[4];
			param_max[4]=param[4];
			cout << "alphas will remain fixed at " << param[1] << " and " << param[4] << endl;
			outfile << "alphas will remain fixed at " << param[1] << " and " << param[4] << endl;
			break;
		case 2:
			param_min[2]=param[2];
			param_max[2]=param[2];
			param_min[5]=param[5];
			param_max[5]=param[5];
			cout << "betas will remain fixed at " << param[2] << " and " << param[5] << endl;
			outfile << "betas will remain fixed at " << param[2] << " and " << param[5] << endl;
			break;
		case 3:
			param_min[1]=param[1];
			param_max[1]=param[1];
			param_min[4]=param[4];
			param_max[4]=param[4];
			cout << "alphas will remain fixed at " << param[1] << " and " << param[4] << endl;
			outfile << "alphas will remain fixed at " << param[1] << " and " << param[4] << endl;
			param_min[2]=param[2];
			param_max[2]=param[2];
			param_min[5]=param[5];
			param_max[5]=param[5];
			cout << "betas will remain fixed at " << param[2] << " and " << param[5] << endl;
			outfile << "betas will remain fixed at " << param[2] << " and " << param[5] << endl;
			break;
		case 4:
			param_min[1]=param[1];
			param_max[1]=param[1];
			param_min[2]=param[2];
			param_max[2]=param[2];
			cout << "alpha0 will remain fixed at " << param[1] << "; beta0 will remain fixed at " << param[2] << endl;
			outfile << "alpha0 fixed at " << param[1] << "; beta0 fixed at " << param[2] << endl;
			break;
		default:
			cout << "invalid selection - both alpha and beta will vary" << endl;
			outfile << "invalid selection - both alpha and beta will vary" << endl;
	}

// do{ //ignore the reduced chi^2 goal for now
// do initial aod calculations to get initial reduced chi^2 values
// start with a random swarm
do{ // run multiple times (populations)
	cout << "initializing swarm..." << endl;
	init_swarm(swarm,param_min,param_max,param,particles,param_fix); // initialize the swarm
	cout << "initializing aods..." << endl;
	for(i=0;i<particles;i++){ // cycle through particles
		for(j=0;j<channels;j++){ // cycle through wavelengths
			calc_aod(swarm[i],mie[j],aodc[i][j],dist); // calculate initial aods with first random swarm
		}
		switch(param_fix)	// instead of setting param_min and param_max to be equal to param and redoing calculations, just skip updates for fixed parameters
		{
/*			case 0:
				for(j=0;j<6;j++){
					Pl[i][j]=swarm[i][j]; // initialize local minima for N0, a, b using first random swarm
				}
				break; */
			case 1:
				Pl[i][0]=swarm[i][0]; // initialize local minima for N0 using first random swarm
				Pl[i][1]=param[1]; // initialize local minima for a using fixed value
				Pl[i][2]=swarm[i][2]; // initialize local minima for b using first random swarm
				Pl[i][3]=swarm[i][3]; // initialize local minima for N1 using first random swarm
				Pl[i][4]=param[4]; // initialize local minima for a1 using fixed value
				Pl[i][5]=swarm[i][5]; // initialize local minima for b1 using first random swarm
				break;
			case 2:
				Pl[i][0]=swarm[i][0]; // initialize local minima for N0 using first random swarm
				Pl[i][1]=swarm[i][1]; // initialize local minima for a using first random swarm
				Pl[i][2]=param[2]; // initialize local minima for b using fixed value
				Pl[i][3]=swarm[i][3]; // initialize local minima for N1 using first random swarm
				Pl[i][4]=swarm[i][4]; // initialize local minima for a1 using first random swarm
				Pl[i][5]=param[5]; // initialize local minima for b1 using fixed value
				break;
			case 3:
				Pl[i][0]=swarm[i][0]; // initialize local minima for N0 using first random swarm
				Pl[i][1]=param[1]; // initialize local minima for a using fixed value
				Pl[i][2]=param[2]; // initialize local minima for b using fixed value
				Pl[i][3]=swarm[i][3]; // initialize local minima for N1 using first random swarm
				Pl[i][4]=param[4]; // initialize local minima for a1 using fixed value
				Pl[i][5]=param[5]; // initialize local minima for b1 using fixed value
				break;
			default:
				for(j=0;j<6;j++){
					Pl[i][j]=swarm[i][j]; // initialize local minima for N0, a, b using first random swarm
				}
				break;
		}
	}
/*
//testing loop
	for(i=0;i<channels;i++){
		cout << "aodc[0]["<<i<<"] = " << aodc[0][i] << " aodc[1]["<<i<<"] = " << aodc[1][i] << " aodc[2]["<<i<<"] = " << aodc[2][i] << endl;
	}
*/
	cout << "calculating initial reduced chi^2 values..." << endl;

	calc_red_chi2(aodm, aodc, f0, dist); // find initial reduced chi^2 values
/*
//testing loop
	for(i=0;i<particles;i++){
		cout << "f0["<<i<<"] = " << f0[i] << endl;
	}
*/

	cout << "initializing swarm..." << endl;
	init_swarm(swarm,param_min,param_max,param,particles,param_fix); // initialize the swarm
/*	for(i=0;i<particles;i++){
		cout << "Particle " << i << ": " << swarm[i][0] << ", " << swarm[i][1] << ", " << swarm[i][2] << endl;
	}
*/
	cout << "Population \t Generation \t N_t" << endl;
	while(gen<gen_max){ // let the swarm evolve
		cout << population << "\t" << gen << "\t" << N_t << endl;
//		cout << "calculating aods..." << endl;
		for(i=0;i<particles;i++){ // cycle through particles
			for(j=0;j<channels;j++){ // cycle through wavelengths
				calc_aod(swarm[i],mie[j],aodc[i][j],dist); // aods calculated with individual particles' N0, a, b
			}
		}
//		cout << "calculating reduced chi^2..." << endl;
		calc_red_chi2(aodm, aodc, f, dist); // find reduced chi^2 values
//		cout << "finding local minima..." << endl;
		find_local_min(f, f0, Pl, swarm); // find local minima
//		cout << "finding global minimum..." << endl;
		i_g = find_global_min(f, Pg, swarm, particles); // find global minimum
//		cout << "updating global minimu..." << endl;
		for(i=0;i<6;i++){ // cycle through parameters N0, a, b and
			Pg[i]=swarm[i_g][i]; // set global minima for N0, a, b
		}
		init_swarm(Pj,param_min,param_max,param,scouts,param_fix); // generate Pj particle for improving global minimum searching
		for(i=0;i<scouts;i++){ // cycle through each scout
			for(j=0;j<channels;j++){ // cycle through wavelengths
				calc_aod(Pj[i], mie[j], aodj[i][j],dist); // aods for extra Pj particles
			}
		}
		calc_red_chi2(aodm,aodj,fj, dist); // reduced chi^2 for extra Pj particle
		i_j = find_global_min(fj,Pj0,Pj,scouts);
		if(fj[i_j]<f[i_g]){ // if extra Pj particles' best reduced chi^2 is less than current global minimum, make it the new global minimum
			for(i=0;i<6;i++){ // cycle through parameters N0, a, b
				Pg[i]=Pj0[i];
			}
			cout << "Scouts helped!" << endl;
			outfile << "Scouts helped! Generation: " << gen << endl;
		}
//		cout << "updating the swarm..." << endl;
		update_swarm(swarm, Pl, Pg, param_min, param_max, param, param_fix);
//		cout << "checking stop condition..." << endl;
		d_global_small=abs(Pg[0]-Pg0[0])/Pg0[0] < 1.0E-5 && abs(Pg[1]-Pg0[1])/Pg0[1] < 1.0E-5 && abs(Pg[2]-Pg0[2])/Pg0[2] < 1.0E-5;
		if(dist==4){d_global_small=d_global_small && abs(Pg[3]-Pg0[3])/Pg0[3] < 1.0E-5 && abs(Pg[4]-Pg0[4])/Pg0[4] < 1.0E-5 && abs(Pg[5]-Pg0[5])/Pg0[5] < 1.0E-5;}
		if (d_global_small){
			N_t+=1; // if no significant change update number of trials for no change
		}
		else {
			N_t=0; // if significant change reset number of trials
		}
		if(N_t>N_tmax){ // stop condition reached?
			cout << "delta N0: " << (Pg[0]-Pg0[0])/Pg0[0] << "  delta a: " << (Pg[1]-Pg0[1])/Pg0[1] << "  delta b: " << (Pg[2]-Pg0[2])/Pg0[2] << endl;
			if(dist==4){
				cout << "delta N1: " << (Pg[3]-Pg0[3])/Pg0[3] << "  delta a1: " << (Pg[4]-Pg0[4])/Pg0[4] << "  delta b1: " << (Pg[5]-Pg0[5])/Pg0[5] << endl;
			}
			cout << "reduced chi^2 = " << f[i_g] << endl;
			break; // if stop condition met, break out of evolution loop			
		}
		else{
			for(i=0;i<6;i++){ // cycle through parameters N0, a, b
				Pg0[i]=Pg[i]; // to update global min
			}
		}
		gen++;
	} // end of evolution loop
// N_t=0; // reset number of trials
// population += 1;
// gen=0; // reset generation
// } while(f[i_g]>reduced chi^2_goal && population < 11); // limit to 10 populations or desired reduced chi^2 goal

	outfile << Pg[0] << ", " << Pg[1] << ", " << Pg[2] << endl;
	if(dist==4){
		outfile << Pg[3] << ", " << Pg[4] << ", " << Pg[5] << endl;
	}
	outfile << "Generations: " << gen << "\t" << "Population: " << population << endl;
	cout << "Generation # " << gen << ", N0 = " << Pg[0] << ", alpha = " << Pg[1] << ", beta = " << Pg[2] << endl;
	if(dist==4){
		cout << "N1 = " << Pg[3] << ", alpha1 = " << Pg[4] << ", beta1 = " << Pg[5] << endl;
	}
	for(i=0;i<channels;i++){
		cout << "Measured: " << aodm[i] << "  Calculated: " << aodc[i_g][i] << endl;
		outfile << "Measured: " << aodm[i] << "  Calculated: " << aodc[i_g][i] << endl;
	}
	outfile << "reduced chi^2 = " << f[i_g] << endl;
	outfile1 << f[i_g] << ", " << Pg[0] << ", " << Pg[1] << ", " << Pg[2];
	if(dist==4){
		outfile1 << ", " << Pg[3] << ", " << Pg[4] << ", " << Pg[5];
	}
	outfile1 << endl;
	population += 1;
	gen = 0;
	N_t = 0;
} while(population<pop_max+1); // limit to pop_max populations
	outfile.close();
	outfile1.close();
	return 0;
}

double random(double low, double high)
{
	double result,range;
	range = log10(high/low); // log(high)-log(low)=log(high/low)
	result = log10(low) + (rand() % 1000000000)/1.0E9*range;
	return pow(10,result);
}

void calc_aod(vector<double> &part, vector<double> &Qext, double &tauc, int dist)
{
	int i=0;
	double tau1=0.0,tau2=0.0;
	tauc=0;
	while(i<9962){ // cycle through all radii
		if(dist==1){ // Junge distribution
			tauc += pow((i+40)*CONV_FACT2,2.0)*Qext[i]*pow((i+40.5)*CONV_FACT2,-(part[1]+1)); // Junge distribution
		}
		else if(dist==2){ // Gamma distribution
			tauc += pow((i+40)*CONV_FACT2,2.0)*Qext[i]*pow((i+40.5)*CONV_FACT2,part[2])*exp(-part[1]*(i+40.5)*CONV_FACT1); // Gamma distribution
		}
		else if(dist==3){ // Lognormal distribution
			tauc += pow((i+40),2.0)/(i+40.5)*CONV_FACT2*Qext[i]*exp(-pow(log((i+40.5)*CONV_FACT1/part[1])/part[2],2.0)/2.0); // Lognormal distribution
		}
		else if(dist==4){ // Bimodal Lognormal distribution
			tau1 += pow((i+40),2.0)/(i+40.5)*CONV_FACT2*Qext[i]*exp(-pow(log((i+40.5)*CONV_FACT1/part[1])/part[2],2.0)/2.0); // Aitkin mode
			tau2 += pow((i+40),2.0)/(i+40.5)*CONV_FACT2*Qext[i]*exp(-pow(log((i+40.5)*CONV_FACT1/part[4])/part[5],2.0)/2.0); // accumulation mode
		}
		i++;
	}
	if(dist==3){	// lognormal distribution
		tauc *= sqrt(PI/2.0)*1.0E-7*part[0]/part[2]; // multiply by the constants from the sum: sqrt(pi/2), 2*delta r, N0/beta
	}
	else if(dist==4){	// bimodal lognormal distribution
		tauc = tau1*sqrt(PI/2.0)*1.0E-7*part[0]/part[2]+tau2*sqrt(PI/2.0)*1.0E-7*part[3]/part[5];
	}
	else if(dist==2){	// modified gamma distribution
		tauc *= 2.0*PI*5.0E-8*part[0]; // multiply by the constants from the sum: pi, 2*delta r, N0, beta
	}
	else{
		tauc *= 2.0*PI*5.0E-8*part[0]*part[2]*1.0E-12; // multiply by the constants from the sum: pi, 2*delta r, N0, beta and convert from number at 1cm to number at 1um
	}
	return;
}

void init_mie(matrix &mie)
// fills the matrix mie[wavelength][radius]
{
	string fname;
	fname="mie-data/mie412nm.csv";
	initQ(mie[0],fname);
	fname="mie-data/mie440nm.csv";
	initQ(mie[1],fname);
	fname="mie-data/mie463nm.csv";
	initQ(mie[2],fname);
	fname="mie-data/mie479nm.csv";
	initQ(mie[3],fname);
	fname="mie-data/mie500nm.csv";
	initQ(mie[4],fname);
	fname="mie-data/mie520nm.csv";
	initQ(mie[5],fname);
	fname="mie-data/mie556nm.csv";
	initQ(mie[6],fname);
	fname="mie-data/mie610nm.csv";
	initQ(mie[7],fname);
	fname="mie-data/mie675nm.csv";
	initQ(mie[8],fname);
	fname="mie-data/mie750nm.csv";
	initQ(mie[9],fname);
	fname="mie-data/mie778nm.csv";
	initQ(mie[10],fname);
	fname="mie-data/mie870nm.csv";
	initQ(mie[11],fname);
	fname="mie-data/mie1020nm.csv";
	initQ(mie[12],fname);
	return;
}

void initQ(vector<double> &Qext, string fname)
// fills an individual row (i.e., wavelength) of mie[wavelength][radius]
{
	int i;
	double r, Qsca, Qabs, Gsca;	// data in Mie files but not needed
	char *fname1,dataline[100];
	string mystr;
	fstream miefile;
	fname1 = new char [fname.size()+1];	/* c_string method */
	strcpy(fname1,fname.c_str());

	miefile.open(fname1);
	if (miefile.is_open())
	{
		i=0;
		getline(miefile,mystr); // throw away header line
		while(i<9962)
		{
			getline(miefile,mystr);
			strcpy(dataline,mystr.c_str());
			sscanf(dataline,"%lf,%lf,%lf,%lf,%lf",&r,&Qext[i],&Qsca,&Qabs,&Gsca); // fills the Qext[] array ignoring other values
			i++;
		}
		miefile.close();
		cout << "Processed mie file " << fname << "\n";
	}
	else cout << "\nUnable to open file " << fname << "\n";

	return;
}

void calc_red_chi2(double aodm[], matrix &aodc, double f[], int dist)
{
	int i=0,j=0,rows;
	double a;
	rows=aodc.size();
	for(i=0;i<rows;i++){ // cycle through all particles
		a=0.0;
		for(j=0;j<channels;j++){ // cycle through all wavelengths
			a += pow(aodm[j]-aodc[i][j],2.0)/aodc[i][j]; // sum of square deviations
		}
		switch (dist) {
		case 1:		// Junge distribution
			f[i]=a/(channels-2);
			break;
		case 4:		// bimodal lognormal distribution
			f[i]=a/(channels-6);
			break;
		default:	// modified gamma and lognormal distributions
			f[i]=a/(channels-3);
			break;
		}
	}
	return;
}

void init_swarm(matrix &swarm, double param_min[], double param_max[], vector<double> &param, int num, int param_fix)
{
	int i=0;
	while(i<num){ // cycle through all particles
		init_particle(swarm[i], param_min, param_max, param, param_fix); // generate parameters
/*
		cout << "Particle " << i << " initialized: " << swarm[i][0] << ", " << swarm[i][1] << ", " << swarm[i][2] << endl;
*/
		i++;
	}
	return;
}

void init_particle(vector<double> &part, double param_min[], double param_max[], vector<double> &param, int param_fix)
{
	int j;
	switch (param_fix)
	{
/*		case 0: // all parameters vary
			for(j=0;j<6;j++){
				part[j]=random(param_min[j],param_max[j]); // set N0, a, and b parameters
			}
			break; */
		case 1: // alphas are fixed
			part[0]=random(param_min[0],param_max[0]); // set N0 parameter
			part[1]=param[1];	// fix a parameter
			part[2]=random(param_min[2],param_max[2]); // set b parameter
			part[3]=random(param_min[3],param_max[3]); // set N1 parameter
			part[4]=param[4];	// fix a1 parameter
			part[5]=random(param_min[5],param_max[5]); // set b1 parameter
			break;
		case 2: // betas are fixed
			part[0]=random(param_min[0],param_max[0]); // set N0 parameter
			part[1]=random(param_min[1],param_max[1]); // set a parameter
			part[2]=param[2];	// fix b parameter
			part[3]=random(param_min[3],param_max[3]); // set N1 parameter
			part[4]=random(param_min[4],param_max[4]); // set a1 parameter
			part[5]=param[5];	// fix b1 parameter
			break;
		case 3: // alphas and betas are fixed
			part[0]=random(param_min[0],param_max[0]); // set N0 parameter
			part[1]=param[1];	// fix a parameter
			part[2]=param[2];	// fix b parameter
			part[3]=random(param_min[3],param_max[3]); // set N1 parameter
			part[4]=param[4];	// fix a1 parameter
			part[5]=param[5];	// fix b1 parameter
			break;
		case 4: // alpha0 and beta0 are fixed
			part[0]=random(param_min[0],param_max[0]); // set N0 parameter
			part[1]=param[1];	// fix a parameter
			part[2]=param[2];	// fix b parameter
			part[3]=random(param_min[3],param_max[3]); // set N1 parameter
			part[4]=random(param_min[4],param_max[4]); // set a1 parameter
			part[5]=random(param_min[5],param_max[5]); // set b1 parameter
			break;
		default:
			for(j=0;j<6;j++){
				part[j]=random(param_min[j],param_max[j]); // set N0, a, and b parameters
			}
			break;
	}
	return;
}

int find_global_min(double f[], double Pg[], matrix &swarm, int num)
{
	int i;
	int g=0;
	for(i=0;i<num;i++){ // cycle through all particles
		if (f[i]<f[g]) { // look for lowest reduced chi^2 value
			g = i; // set i_g to index of lowest reduced chi^2 value
		}
	}
//	for(i=0;i<3;i++){ // cycle through N0, a, b
//		Pg[i]=swarm[i_g][i]; // reset the global minimum
//	}
	return g; // return index of lowest reduced chi^2 value
}

void find_local_min(double f[], double f0[], matrix &Pl, matrix &swarm)
{
	int i,j;
	for(i=0;i<particles;i++){ // cycle through all particles
		if(f[i]<f0[i]){ // test if current reduced chi^2 is less than previous lowest reduced chi^2
			f0[i]=f[i]; // if so, reset the lowest reduced chi^2 and
			for(j=0;j<6;j++){ // reset the local minimum
				Pl[i][j]=swarm[i][j];
			}
			j=0;
		}
	}
	return;
}

void update_swarm(matrix &swarm, matrix &Pl, double Pg[], double param_min[], double param_max[], vector<double> &param, int param_fix)
/* 
Implements Eq (10) from Yuan's article:
X_i(t+1)=X_i(t)+c1*r1*[P_i(t)-X_i(t)]+c2*r2*[P_g(t)-X_i(t)]
*/
{
	int i,j; // particle and parameter indices
	double r1,r2; // random numbers in interval [0,1]
	bool out_of_bounds; // variable to test if particle is out of the domain
	i=0;
	while(i<particles){
		out_of_bounds=false;
		r1 = rand() % 1000000000/1.0E9; // each particle gets a new random number
		r2 = rand() % 1000000000/1.0E9;
		switch (param_fix)
		{
			case 1: // only update N and beta
				swarm[i][0] = (1.0-c1*r1-c2*r2)*swarm[i][0]+c1*r1*Pl[i][0]+c2*r2*Pg[0];
				swarm[i][2] = (1.0-c1*r1-c2*r2)*swarm[i][2]+c1*r1*Pl[i][2]+c2*r2*Pg[2];
				swarm[i][3] = (1.0-c1*r1-c2*r2)*swarm[i][3]+c1*r1*Pl[i][3]+c2*r2*Pg[3];
				swarm[i][5] = (1.0-c1*r1-c2*r2)*swarm[i][5]+c1*r1*Pl[i][5]+c2*r2*Pg[5];
				break;
			case 2: // only update N and alpha
				swarm[i][0] = (1.0-c1*r1-c2*r2)*swarm[i][0]+c1*r1*Pl[i][0]+c2*r2*Pg[0];
				swarm[i][1] = (1.0-c1*r1-c2*r2)*swarm[i][1]+c1*r1*Pl[i][1]+c2*r2*Pg[1];
				swarm[i][3] = (1.0-c1*r1-c2*r2)*swarm[i][3]+c1*r1*Pl[i][3]+c2*r2*Pg[3];
				swarm[i][4] = (1.0-c1*r1-c2*r2)*swarm[i][4]+c1*r1*Pl[i][4]+c2*r2*Pg[4];
				break;
			case 3: // only update N
				swarm[i][0] = (1.0-c1*r1-c2*r2)*swarm[i][0]+c1*r1*Pl[i][0]+c2*r2*Pg[0];
				swarm[i][3] = (1.0-c1*r1-c2*r2)*swarm[i][3]+c1*r1*Pl[i][3]+c2*r2*Pg[3];
				break;
			case 4: // only update N0,N1,a1,b1
				swarm[i][0] = (1.0-c1*r1-c2*r2)*swarm[i][0]+c1*r1*Pl[i][0]+c2*r2*Pg[0];
				swarm[i][3] = (1.0-c1*r1-c2*r2)*swarm[i][3]+c1*r1*Pl[i][3]+c2*r2*Pg[3];
				swarm[i][4] = (1.0-c1*r1-c2*r2)*swarm[i][4]+c1*r1*Pl[i][4]+c2*r2*Pg[4];
				swarm[i][5] = (1.0-c1*r1-c2*r2)*swarm[i][5]+c1*r1*Pl[i][5]+c2*r2*Pg[5];
				break;
			default: // update all
				j=0;
				while(j<6){
					swarm[i][j] = (1.0-c1*r1-c2*r2)*swarm[i][j]+c1*r1*Pl[i][j]+c2*r2*Pg[j];
					// should this only generate a new particle after j=5?
/*					if (swarm[i][j] < param_min[j] || swarm[i][j] > param_max[j]){ // check to see if out bounds
						init_particle(swarm[i], param_min, param_max, param_fix); //if so, generate random postion in bounds
					} */
					j++;
				}
				break;
		}
		for(j=0;j<6;j++){
			out_of_bounds=out_of_bounds || (swarm[i][j]<param_min[j]) || (swarm[i][j]>param_max[j]); // test for parameters out of bounds
		}
		if(out_of_bounds){init_particle(swarm[i], param_min, param_max, param, param_fix);} // if any parameter out of bounds, generate new particle
		i++;
	}
	return;
}