// each call to forward produces signals for all nch channels
#include "FORTRANarray.hpp"
#include "Level1Data.h"
#include "fortranfunctions.h"
#include "GATS_Utilities.hpp"

#include <iostream>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <cassert>
#include <fstream>

using namespace GATS_Utilities;

float CalcDiff(float p, const float deltap, const float pmin,  const float pmax)
{
	p+= deltap;
	if( p<= pmax) return p;
	std::cout << p << "  " << deltap << "  " << pmin << std::endl;
	p = 0.5*p;
	std::cout << p << std::endl;
	assert( p>= pmin);
	return p;
}


extern "C" void derivative( 
	int* ilev,	int* igdo, 		int* ngdo,
	int* igas,	int* ngas,		int* num_iterations,
	int* nl,	int* nch,		int* itab,
	int *nreg,	int ***ievent,		int* maxit,
	float ***vapr,	float **vaprt,		float* tst,
	float ***varm,	double **meas_signals,
	float **evar,	float* zt,
	float* pa,	float* ta,	float* za,
	float **qmix,	double **sim_signals,	
	int* ihydro,	int* first_call,	int* converged,
	int* bottom_up, int* nrays, int *iscale	 )
{
	const static float FMIN = 1.e-25;
	const static double DMAX = 1.e25;
	const static double DMIN = 1.e-25;
	const static double TESTF = 3.0 ;  // derivative must cause at least a 3 noise unit change.
//	const static double TSTF = 0.01;
	const static int MAXDO = 9;
	const static float CHG = 0.02;
	const static float CHGT = 1.0;
	const static float T1 = 110.0;
	const static float T2 = 990.0;
	const static float Q1 = 1.e-14;
	const static float Q2 = 1.0;
	

	static FORTRANarray<float> Qold ,Q_apriori;
	static FORTRANarray<double> MeasSN, Vary, Wght ;
 	static std::vector<float> Pold, Told,  P_apriori, T_apriori;
	static std::vector<int> igasIndex, iGo, iDerivCount;
	static std::vector<double> DeltaQ,  SimSN, SigOld;
	static bool NewUpdate;
	static bool RetrievalIsGone;
	static bool DerivativeIsJunk=false;
	static bool scaleAbove;
	static int ido;
	static int l_nrays; 

	float fact;
	int ik;
	
	
	//std::cout << " nrays " << *nrays << "  nl " << *nl << std::endl;  exit(23);
/*
static bool once =true;
static std::ofstream dfile;

if(once) {
  
 	dfile.open("derivatives.txt");

	for(int il = 0; il< *nl; ++il) {
		dfile << ta[il] << "   " << za[il] << std::endl;
		//dfile << qmix[0][il] << "  " << zt[il] <<  std::endl;
	}
	//dfile << "&" << std::endl;
	once=false;
}
*/ 

// initial call for the current REtrieval setup


/*
     static int badass= 1;
     *converged = 0;

 std::cout << badass << "  " << 
    ConvertToStringPrec(meas_signals[0][*ilev-1]) << "  " <<
    ConvertToStringPrec(sim_signals[0][*ilev-1])  << "  " << za[*ilev-1] << std::endl ;

     qmix[0][*ilev-1] = 1.e-6;

     if(badass %3 ==0 ) qmix[0][*ilev-1] =  2.e-6; 
     if(badass ==8 ) exit(23);
     badass++;
     return;
*/     
////////////////////////////////////////////////////////////////



	if(*first_call != 0) {
		assert( *nch >= 1  );  // need at least one channel
		assert( *nl >= 1 ); 
		assert( *ngdo >= 1 );

		assert( *ngdo <= *ngas+1 );  
// allocate static arrays	
		iGo = std::vector<int>( *nch);
		iDerivCount = std::vector<int>( *ngdo,0);
		DeltaQ = std::vector<double>(*ngdo);
		Wght = FORTRANarray<double>(*ngdo, *nch);
		MeasSN = FORTRANarray<double> ( *nl, *nch);

		Vary = MeasSN;
//std::cout << "********************** " << std::endl;
		for(int il=0; il < *nl; ++il) {
			int idet = itab[il]-1;
			for(int ich=0; ich < *nch; ++ich) {
				if(varm[ich][idet][il] > FMIN) {
					Vary.setVal( 1.0/  varm[ich][idet][il] , il+1, ich+1);
				} else {
					Vary.setVal( DMAX, il+1, ich+1 );
				}
				MeasSN.setVal( meas_signals[ich][il] * std::sqrt( Vary(il+1,ich+1) ), il+1, ich+1);

//std::cout << 1.0/std::sqrt( Vary(il+1,ich+1) )  << "  " << za[il] << "  " << std::endl; //meas_signals[ich][il] << "  " << MeasSN(il+1, ich+1)  << std::endl;
//std::cout << za[il] << "  " << meas_signals[ich][il] << "  " << MeasSN(il+1, ich+1)  << std::endl;
			}
		}
// exit(23);
//std::cout << "********************** " << std::endl;

		SimSN = std::vector<double>(*nch);
		SigOld = SimSN;
	
		Qold = FORTRANarray<float>( *nl, *ngas);
		Pold = std::vector<float>( *nl);
		Told = std::vector<float>( *nl);
		std::copy( ta, ta+ *nl, Told.begin() );	
		std::copy( pa, pa+ *nl, Pold.begin() );
		for(int ig=0; ig < *ngas; ++ig) {
std::cout << "FIRST CALL  igas " << igas[ig] << std::endl;
			std::copy(&qmix[ig][0], &qmix[ig][0]+ *nl, Qold.getArray(1,ig+1) );
		}


//     for (int i = 0; i< *nl; ++i ) 
//      std::cout <<  ta[i] << "   " << za[i] << std::endl;
//endl   std::cout << "&" << std::endl;
/*
for (int j=0; j< *ngas; ++j) {
std::cout << "gas " <<  igas[j] << std::endl;
for (int i = 0; i< *nl; ++i ) {
std::cout << qmix[j][i] << " " << Qold(i+1,j+1) << "  " << za[i] << std::endl;	
}
}
exit(23) ;
*/		
		T_apriori = Told;
		P_apriori = Pold;
		Q_apriori = Qold;

		igasIndex.clear();
		for(int i = 0; i< *ngdo; ++i) {
			assert( igdo[ i] != 0);   // what's a zero doing here?
			if( igdo[ i] > 0 ) {
				int* p = std::find( igas, igas+ *ngas, igdo[i] );
				int k = std::distance(igas,p);
				assert(k < *ngas);
				igasIndex.push_back(k);
			} else {
				igasIndex.push_back( -1 );   // indicates Temperature
			}

//			std::fill(&evar[i][0], &evar[i][*nl], -999.0);
            std::fill(&evar[i][0], &evar[i][*nl], 0.0);

		}

		l_nrays = (*bottom_up ) ? *nrays : 1;

//std::cout << *bottom_up << "  " << *iscale << "  " << *nrays << std::endl; exit(23);
		if(! *bottom_up && *iscale) {
			l_nrays=1000;
			scaleAbove=true;
			DerivativeIsJunk=false;
		} else {
			scaleAbove=false;
		}
		assert( l_nrays > 0);
		ido = 0;
		NewUpdate= true;
		RetrievalIsGone = false;
		*first_call = 0;
 std::cout << " ok finished first call " << std::endl;
	}  //completed with first call


	*converged = 0;   //set convergence to 0  or false
	*ihydro = 0;   // don't do hydrostatics unless we have to (determined below)

	if(NewUpdate) {
		for(int ich=0; ich < *nch; ++ich) {
			SimSN[ich] = sim_signals[ich][*ilev-1] * std::sqrt( Vary( *ilev, ich+1) );
//            std::cout << *ilev << " " << sim_signals[ich][*ilev-1]  
//                      << "    " <<  meas_signals[ich][*ilev-1] << std::endl;
		}
                SigOld = SimSN;

//std::cout << "In update " << ConvertToStringPrec(MeasSN(*ilev, 1)) << "   " << ConvertToStringPrec(SimSN[0] ) << std::endl;

//std::cout << " in update signal is " << ConvertToStringPrec(SigOld[0]) << "  " << za[*ilev-1] << std::endl;

		std::copy( ta, ta+ *nl, Told.begin() );	
		std::copy( pa, pa+ *nl, Pold.begin() );

		for(int ig=0; ig < *ngas; ++ig) {
			std::copy(&qmix[ig][0], &qmix[ig][0]+ *nl, Qold.getArray(1,ig+1) );
		}
		
std::cout << " we update " << qmix[0][*ilev-1 ] << "  " << *num_iterations << std::endl;
//std::cout << " Ok we are in Update " << *num_iterations <<  "  " << za[*ilev-1] << std::endl;
		NewUpdate= false;
	}
//std::cout << " this is ngdo " << *ngdo << std::endl;

	if( ido == *ngdo ) goto skip_derivative;
	
std::cout << " do derivative " << qmix[0][*ilev-1] << std::endl;
// Derivative calculation section
derivative:
	ik = igasIndex[ido];
	float factr;
	if(ik >=0 ) {
       	int idg = (igas[ik] > 10000) ? igas[ik]/10000 : igas[ik];
		factr = ( idg < 9000) ? 4.5 : 6.0;   // idg >=9000 indicates aerosol extinction
	} else {
		factr = 2.0;
	}
// We have no derivative attempts for this ido so we need to compute one
	if( iDerivCount[ido] == 0) {   
		iDerivCount[ido]++;
 // If retrieving temperature.. never change the temperature at the reference altitude.
		if( (ik < 0 &&  *ilev == *nreg)  || RetrievalIsGone   ) {  
			std::fill(iGo.begin() , iGo.end() , 0);
			for(int ich=0; ich< *nch; ++ich) Wght.setVal(0.0,ido+1,ich+1);
		} else {
			std::fill(iGo.begin() , iGo.end() , 1);
			if( ik >= 0) {
				qmix[ik][*ilev -1] = 
					CalcDiff( Qold( *ilev, ik+1), CHG * Qold( *ilev, ik+1) , Q1,  Q2);
				fact = qmix[ik][*ilev -1] / std::max(Q1,Qold( *ilev, ik+1));
std::cout << Qold( *ilev, ik+1) << "  " << CHG * Qold( *ilev, ik+1) << "  " << qmix[ik][*ilev -1] << std::endl; exit(23);
//std::cout << " lrays " << l_nrays << " fact " << fact << std::endl; exit(23);

				for(int il= *ilev-1; il >= std::max(1, *ilev -l_nrays + 1) ; --il ) {
					qmix[ik][il -1] = fact * Qold(il, ik+1);
					if(qmix[ik][il-1] < Q1) qmix[ik][il-1] = Q1;
					if(qmix[ik][il-1] > Q2) qmix[ik][il-1] = Q2;
				}
			} else {  //this is a temperature change so we call hydrostatics
				ta[ *ilev-1] = 
					CalcDiff(  Told[ *ilev-1],   CHGT , T1,  T2);
				fact = ta[*ilev-1] - Told[*ilev-1];
				for(int il= *ilev-1; il >= std::max(1, *ilev -l_nrays + 1) ; --il ) {
					ta[il-1] = Told[ il-1] + fact;
					if(ta[il-1] < T1) ta[il-1] = T1;
					if(ta[il-1] > T2) ta[il-1] = T2;
				}
				*ihydro = 1;
			}
                   std::cout << qmix[ik][*ilev-1] << "  " << Qold( *ilev, ik+1) << std::endl; exit(23);
			return;  // return to do a forward with this update applied 
		}
// We gave ido a perturbation now we see if that led to a good derivative.
	} else {
//std::cout << " first iteration " << std::endl;
		for(int ich=0; ich< *nch; ++ich) {
//std::cout << " in loop " << *nch << std::endl;
			if(ievent[ich][ido][ *ilev-1] != 0 ) {
//std::cout << "in if statement " << std::endl;
				if( iGo[ich] == 0) continue;   // we already have an ido derivative for this channel
				double deltaP;
				if(ik >= 0) {
					deltaP = qmix[ik][*ilev -1] - Qold(*ilev, ik+1);
				} else {
					deltaP = ta[*ilev -1] - Told[*ilev-1];
				}

				double signalNew = sim_signals[ich][*ilev-1] * std::sqrt( Vary( *ilev, ich+1) );
//					double testF = std::fabs(signalNew * TSTF);
				double deltaS = signalNew - SimSN[ich];

//std::cout << za[*ilev-1] << "  " << sim_signals[ich][*ilev-1] << "  " << meas_signals[ich][*ilev-1] << 
//" " << deltaS <<   "  " << deltaP  << std::endl; 

// We have an adequate change in signal so we create the derivative
				if(std::fabs(deltaS) > TESTF) {
					Wght.setVal( deltaS/deltaP, ido+1, ich+1 );
					iGo[ich] = 0;
					if(ik >=0) {    // we reset the ido profile to its original unperturbed state
						std::copy(Qold.getArray(1,ik+1),Qold.getArray(*nl,ik+1)+1,&qmix[ik][0]);
					} else {
						std::copy( Told.begin(), Told.end(), ta );	
						std::copy( Pold.begin(), Pold.end(), pa );
						*ihydro=1;
					}
					DerivativeIsJunk= false;
// after maxdo attempts to get a derivative with no success, we set the derivative = 0 if we have done no
// solution iterations or keep our last good one if we have at least made a few guesses.
				} else if( iDerivCount[ido] == MAXDO) {
					iGo[ich] = 0;
					if(*num_iterations == 0) Wght.setVal(0.0,ido+1,ich+1);
					if(ik >=0) {  // again we reset the ido profile to the orignal unperturbed state
						std::copy(Qold.getArray(1,ik+1),Qold.getArray(*nl,ik+1)+1,&qmix[ik][0]);
					} else {
						std::copy( Told.begin(), Told.end(), ta );	
						std::copy( Pold.begin(), Pold.end(), pa );
						*ihydro=1;
					}
					if(*bottom_up) {
						RetrievalIsGone = true;
//						std::cout << " retrieval is gone at " << za[*ilev-1] << std::endl;
					} 
                    			DerivativeIsJunk=true;
					
// we apply a different perturbation to try to get a derivative.
				} else  {
					float delta;
					if( ik >= 0) {
						delta = Qold(*ilev, ik+1) * CHG * 
							std::pow(factr,iDerivCount[ido] );
                                                if(delta > 0.1*Q2) delta = 0.1*Q2;
						qmix[ik][*ilev -1] = 
							CalcDiff( Qold( *ilev, ik+1), delta , Q1,  Q2);
						fact = qmix[ik][*ilev -1] / std::max(Q1,Qold( *ilev, ik+1));
						for(int il= *ilev-1; il >= std::max(1, *ilev -l_nrays +1) ; --il ) {
							qmix[ik][il -1] = fact * Qold(il, ik+1);
							if(qmix[ik][il-1] < Q1) qmix[ik][il-1] = Q1;
							if(qmix[ik][il-1] > Q2) qmix[ik][il-1] = Q2;
						}
					} else {
						delta = CHGT * std::pow(factr,iDerivCount[ido] );
                                                if(delta > 0.1*T2) delta = 0.1*T2;
						ta[ *ilev-1] = 
							CalcDiff(  Told[ *ilev-1],  delta , T1,  T2);
						fact = ta[*ilev-1] - Told[*ilev-1];
						for(int il= *ilev-1; il >= std::max(1, *ilev -l_nrays + 1) ; --il ) {
							ta[il-1] = Told[ il-1] + fact;
							if(ta[il-1] < T1) ta[il-1] = T1;
							if(ta[il-1] > T2) ta[il-1] = T2;
						}
						*ihydro = 1;
					}

				}
// this channel / gas combination has no sensitivity so we set its derivative to zero
			} else {   
				iGo[ich] = 0;
				Wght.setVal(0.0,ido+1,ich+1);
			}
		} // for loop over ich
  
// Check that we have derivatives for every channel channel for this ido(gas, extinction, temperature)
		if(std::accumulate(iGo.begin(), iGo.end(),0) > 0 ) {
			iDerivCount[ido]++;
			return;
		}
	}
	ido++;
	if (ido < *ngdo) goto derivative;

    SimSN = SigOld ;
    
//    std::cout << " derivative  done " << qmix[0][*ilev-1] << "  " << Wght(1,1) << 
//    "  " << ConvertToStringPrec(SimSN[0]) << "  " << ConvertToStringPrec(SigOld[0]) << 
//    std::endl;

// Completed the derivatives   Yippee!!!!!!!!!!

skip_derivative:

//if( *num_iterations == 1) std::cout << zt[*ilev-1] << "  " << Wght(1,1) << "   " << *num_iterations << std::endl;

	NewUpdate = true;
	std::fill(iDerivCount.begin(), iDerivCount.end(), 0);

// Compute a priori variances
	std::vector<double> sumDx(*ngdo, 0.0);
	FORTRANarray<double> VarX( *ngdo, *ngdo);
	std::fill(VarX.getArray(), VarX.getArray()+VarX.size(), 0.0);



	for(int ig=0; ig< *ngdo; ++ig) {
		ik = igasIndex[ig];
		float vap, delta;
		if(ik >=0) {
			vap = vapr[ik][ *ilev-1][*ilev-1];
			delta = Q_apriori(*ilev, ik+1) - qmix[ik][*ilev-1];

		} else {
			vap = vaprt[*ilev-1][*ilev-1];
			delta = T_apriori[*ilev-1] - ta[*ilev-1];
		}


		double varxi = (vap > FMIN) ? 1.0/ vap : DMAX;
		
//std::cout << za[*ilev-1] << "   here is delta " << delta << " and varxi " << varxi << std::endl;		
		sumDx[ig] = varxi*delta;
		VarX.setVal( varxi, ig+1, ig+1);
//	VarX.setVal( 1.0, ig+1, ig+1);	
		
//std::cout << za[*ilev-1] << "   delta " << delta << " sumdx " << sumDx[ig] << " varxi " << varxi <<  std::endl;		
	}
	
	
	
//

	for(int ich=0; ich< *nch; ++ich) {
		double delr = MeasSN(*ilev, ich+1) - SimSN[ich];
//std::cout << "iterations  delr " << *num_iterations << "    " << delr << "  " << Wght(1,1) << std::endl;
		for(int ig=0; ig < *ngdo; ++ig) {
			double vector = Wght(ig+1, ich+1);
			sumDx[ig] += ( vector * delr);
//std::cout << " add " << 1.0/( vector * delr) << std::endl;

			for(int ig2=0; ig2 < *ngdo; ++ig2) {
//std::cout <<  " here is vector * wgh " << vector*Wght(ig2+1,ich+1) << std::endl;

				VarX.setVal ( ( VarX(ig+1,ig2+1) +  vector*Wght(ig2+1,ich+1) ), ig+1, ig2+1 ) ;
			}
		}
	}


// invert VarX Matrix
	FORTRANarray<double> VinV(*ngdo, *ngdo);
	if( *ngdo > 1) {
		double det;
		std::vector<int> lindex(*ngdo);
		::minv_( VarX.getArray(), VinV.getArray(), ngdo, &det,  &lindex[0]  );
	} else {
	
//	std::cout << " varx " << VarX(1,1) << " " << DMIN << std::endl;  exit(23);
		if(VarX(1,1) > DMIN) {
			VinV.setVal(1.0/VarX(1,1),1,1 );
		} else {
			VinV.setVal(0.0,1,1);
		}
	}

// Compute the Deltas needed for the next Guess and solution variance

	std::fill(DeltaQ.begin(), DeltaQ.end() , 0.0);
	for(int ig=0; ig< *ngdo; ++ig) {
		for(int ig2=0; ig2< *ngdo; ++ig2) {
			if(ig == ig2 && VinV(ig+1,ig2+1) < DMIN) {
				VinV.setVal(0.0,ig+1,ig2+1);
			} else if( VinV(ig+1, ig2+1) > DMAX) {
				VinV.setVal(DMAX, ig+1, ig2+1);
			}
//std::cout << "where deltas are " << za[*ilev-1] << "  " << VinV(ig+1, ig2+1) << 
 //   "   " <<  sumDx[ig2] << "  " <<     VinV(ig+1, ig2+1) * sumDx[ig2] << std::endl;
    
    
			DeltaQ[ig] += (VinV(ig+1, ig2+1) * sumDx[ig2]);
		}
		
// Solution variance.... needs attention for multiple gases and channels		
		evar[ig][*ilev-1] = (float)VinV(ig+1, ig+1);
//std::cout << " delta q " << DeltaQ[ig] << "  " << std::sqrt( evar[ig][*ilev-1]  ) << std::endl;
	}


//  std::cout << za[*ilev-1] << "  " <<  *num_iterations << " Deltaq " << DeltaQ[0] << "  " <<
//  std::sqrt( evar[0][*ilev-1] ) <<  "   " << tst[0] <<  std::endl;	



// check if another iteration is needed
	std::vector<int> iConvg(*ngdo,1);

	for(int ig=0; ig< *ngdo; ++ig) {
		double testc = std::sqrt( evar[ig][*ilev-1] ) * tst[ig];
//std::cout << " test " << testc << "   " << test[ig] << std::endl;
		if(std::fabs(DeltaQ[ig]) <= testc) {
			iConvg[ig] = 0;
		}
	}

	*num_iterations = *num_iterations + 1;
// we have a solution....

//if(*num_iterations == 7 ) exit(23);
	if(std::accumulate(iConvg.begin(), iConvg.end(), 0) == 0)  {
		*converged = 1;
// we reached max iterations
	} else if(*num_iterations > *maxit) {  
		*converged=1;

/*
if(igasIndex[0] < 0) {
     std::cout << " Max Iterations: layer " << *ilev << " "  << 
	   za[*ilev-1] << " " << Told[*ilev-1] << "last Delta " << DeltaQ[0] << 
 "   derivative dT/dS " << -1.0/Wght(1,1) << "   " <<   MeasSN(*ilev, 1) - SimSN[0]
         << "  " << meas_signals[0][*ilev-1]    <<  std::endl;
//	"   derivative dQ/dS " << Wght(1,1) << std::endl;
} else {

    std::cout << " Max Iterations: layer " << *ilev << " "  <<
           za[*ilev-1] << " " << qmix[0][*ilev-1] << " last Delta " << DeltaQ[0] <<
        "   derivative dQ/dS " << -1.0/Wght(1,1) << "   " //<<   MeasSN(*ilev, 1) - SimSN[0]
       << MeasSN(*ilev, 1) << "   " << SimSN[0] 
         << "  " << meas_signals[0][*ilev-1]    <<  std::endl;
}
*/


// apply the new guess see what happens.
	} else  {
		for(int ig=0; ig< *ngdo; ++ig) {
			ik = igasIndex[ig];

			if(ik >= 0) {
			//	float qmin = (Qold(*ilev, ik+1) ) * 0.5;
			//	float qmax = (Qold(*ilev, ik+1) + 1.0) * 0.5;
//std::cout << " Here is our iteration " << Qold(*ilev, ik+1)  << "  " << DeltaQ[ig]  << "  " << za[*ilev-1] << std::endl;
				qmix[ik][*ilev-1] = Qold(*ilev, ik+1) + DeltaQ[ig] ;
				
				
//		std::cout << " apply guess  " << za[*ilev-1] << "   " << 	qmix[ik][*ilev-1] << std::endl;	
				
				if(qmix[ik][*ilev-1] < Q1) qmix[ik][*ilev-1] = Q1;
				if(qmix[ik][*ilev-1] > Q2) qmix[ik][*ilev-1] = Q2;

//				if(qmix[ik][*ilev-1] < qmin) qmix[ik][*ilev-1] = qmin;
//				if(qmix[ik][*ilev-1] > qmax) qmix[ik][*ilev-1] = qmax;

				fact = qmix[ik][*ilev -1] / std::max(Q1,Qold( *ilev, ik+1));
				
//	std::cout << " here is lrays " << l_nrays << std::endl;			
				for(int il= *ilev-1; il >= std::max(1, *ilev -l_nrays + 1) ; --il ) {
					qmix[ik][il -1] = fact * Qold(il, ik+1);
					if(qmix[ik][il-1] < Q1) qmix[ik][il-1] = Q1;
					if(qmix[ik][il-1] > Q2) qmix[ik][il-1] = Q2;
				}
			} else {
//				float tmin = (Told[*ilev-1] + T1) * 0.5;
//				float tmax = (Told[*ilev-1] + T2) * 0.5;
				ta[*ilev-1] = Told[ *ilev-1] + DeltaQ[ig];
				if(ta[*ilev-1] < T1) ta[*ilev-1] = T1;
				if(ta[*ilev-1] > T2) ta[*ilev-1] = T2;
//				if(ta[*ilev-1] < tmin) ta[*ilev-1] = tmin;
//				if(ta[*ilev-1] > tmax) ta[*ilev-1] = tmax;
				fact = ta[*ilev-1] - Told[*ilev-1];
				for(int il= *ilev-1; il >= std::max(1, *ilev -l_nrays + 1) ; --il ) {
					ta[il-1] = Told[ il-1] + fact;
					if(ta[il-1] < T1) ta[il-1] = T1;
					if(ta[il-1] > T2) ta[il-1] = T2;
				}
				*ihydro = 1;
			}
		}
		if(*num_iterations % 3 == 0) ido = 0;   // recalculate derivatives every third iteration
	}
	
//std::cout << za[*ilev-1] << "  " << *converged << std::endl;	
	if ( *converged == 1 ) {
	
	
//	if (*ilev > 105 ) exit(23)  ;  //{   std::cout << "converged " << exit(23);
//	std::cout << "**converged " << za[*ilev-1] << "   " << qmix[0][*ilev-1] << std::endl;



//std::cout << "za " << za[*ilev-1] << "  " << Wght(1, 1) << std::endl;
//       evar[0][*ilev-1] = Wght(1, 1);
//std::cout << za[*ilev-1] << "  " << -1./Wght(1,1) << std::endl;
		ido = 0;
		if(scaleAbove && !DerivativeIsJunk ) {
			l_nrays=1;
			scaleAbove=false;
		}
	}
}