---

copy-ocsTools.cc

Go to the documentation of this file.

// $Id: ocsTools.cc,v 1.35 2002/06/21 00:52:26 lima Exp $
// #####################################################################
//
//    ocsTools.cc
//
//    Functions from several authors and external libraries. 
//
//    Copyright:
//       Optical Fiber Communications Laboratory (OFCL)
//       Computer Science & Electrical Engineering Department (CSEE)
//       University of Maryland Baltimore County (UMBC)
//
//    The copyright applies to functions written in the OFCL members
// #####################################################################

#include "ocsTools.hh"
#include "ocsConst.hh"

extern ofstream LogFile;


// #################################################################
// #################################################################

// Functions written by Ivan Lima Jr.


//########################################################
//## Get the radius of a circle that contains three points
//#  in the 3-dimensional space (x,y,z), returns the 
//#  central point and, more importantly, returns the 
//#  vector that is normal to the circle whose length
//#  is equal to the total arc length.
//#  It will be used as a penalty extraction to feedback
//#  the PMD compensator.
//## By: Ivan Lima
double GetRadiusCircle(double *PointA, double *PointB, double *PointC)
{
   double vectorAB[4];
   double vectorBC[4];
   double vectorABxBC[4];
   double **matrixA = dmatrix(1,3,1,3);
   double **matrixB = dmatrix(1,3,1,1);   

   for (int ii = 1; ii < 4; ii++) {   
      vectorAB[ii] = PointB[ii] - PointA[ii];
      vectorBC[ii] = PointC[ii] - PointB[ii];      
      vectorABxBC[ii] = vectorBC[ii];
   }   

   CrossProduct(vectorAB, vectorABxBC, 1);

 
   for (int ii = 1; ii < 4; ii++) {
      matrixA[1][ii] = vectorAB[ii];
      matrixA[2][ii] = vectorBC[ii];      
      matrixA[3][ii] = vectorABxBC[ii];   
      matrixB[ii][1] = 0;            
   }   
   
   for (int ii = 1; ii < 4; ii++) {
      matrixB[1][1] += vectorAB[ii]*(PointA[ii]+PointB[ii])/2.;   
      matrixB[2][1] += vectorBC[ii]*(PointB[ii]+PointC[ii])/2.;
      matrixB[3][1] += vectorABxBC[ii]*PointB[ii]; 
   }  
      
   //## Now solve X from A*X = B
   //#  where X is a vector at the center of the circle
   gaussj(matrixA,3,matrixB,1);

   double PointCenter[4];
   for (int ii = 1; ii < 4; ii++) {  
      PointCenter[ii] = matrixB[ii][1];
   }
   
   //## Lines for debugging purpose  
   cout << "DistA_Center "
        << GetDistance2Points(PointA,PointCenter,1) << endl;
   cout << "DistB_Center "
        << GetDistance2Points(PointB,PointCenter,1) << endl;
   cout << "DistC_Center "
        << GetDistance2Points(PointC,PointCenter,1) << endl;        
   //## End of line for debugging purposes
   
   //## The central point will return in PointC
   //## The vector that is normal to the circle 
   //#  is returned in PointA. In order to 
   //#  obtain the polarization dispersion vector 
   //#  from it one needs to divide it by the 
   //#  frequency difference (in Hz) between the original
   //#  PointA and PointC.
   //## The Polarization dispersion vector will 
   //#  return in PointA, but it has to be 
   //#  divided by the frequency difference between
   //#  PointA and PointC. It must be done by the user.
   
   double vectorCenterA[4];
   double vectorCenterB[4];   
   double vectorCenterC[4];   
   for (int ii = 1; ii < 4; ii++) {   
      vectorCenterA[ii] = PointA[ii] - PointCenter[ii];
      vectorCenterB[ii] = PointB[ii] - PointCenter[ii];
      vectorCenterC[ii] = PointC[ii] - PointCenter[ii];            
   }      
   double Radius = GetDistance2Points(PointB,PointCenter,1); 
   double ArcLength;
   ArcLength = Radius*( 
               acos( vectorCenterA[1]*vectorCenterB[1]
                    +vectorCenterA[2]*vectorCenterB[2]
                    +vectorCenterA[3]*vectorCenterB[3] )
              +acos( vectorCenterB[1]*vectorCenterC[1]
                    +vectorCenterB[2]*vectorCenterC[2]
                    +vectorCenterB[3]*vectorCenterC[3] ) );    
   

   //## Set the values to be returned   
   for (int ii = 1; ii < 4; ii++) {  
      PointC[ii] = PointCenter[ii];
      if ( sq(vectorABxBC[1])
          +sq(vectorABxBC[2])
          +sq(vectorABxBC[3]) != 0 ) { 
         PointA[ii] = vectorABxBC[ii]
                     /sqrt( sq(vectorABxBC[1])
                           +sq(vectorABxBC[2])
                           +sq(vectorABxBC[3]) )
                     *ArcLength/(2*pi);      
      }               
      else {
        PointA[ii] = 0;
      }          
   }      
   

        
   return Radius;     
}

//####  Returns the distance between 2 points in 3 dimensions
//## By: Ivan Lima
double GetDistance2Points(double *PointA, double *PointB, int FirstCoord)
{
   double Distance = 0;
   for (int ii = 0; ii < 3; ii++) {     
      Distance += sq(PointA[ii+FirstCoord] - PointB[ii+FirstCoord]);
   }
   return sqrt(Distance);   
}


//#### This function takes care of the periodicity of the value ii
//##   The value of ii must be between (-Periodicity) and (2*Periodicity-1)
//##   The result is a value between 0 and (Periodicity-1)
//## By: Aurenice Lima
int Periodicity(int ii, int Periodicity)
{
   int nn;
   if (ii < 0)
      nn = Periodicity + ii;
   else if (ii >= Periodicity)
      nn = ii - Periodicity;   
   else
      nn = ii;         
   return nn;  
}


//## It produces the angles and the matrices that rotates the 
//#  vector to (-1,0,0)
//## By: Ivan Lima
void TransRzRyV3to_Xhat( double *Vector, double *thetaZthetaY)
//                        ,double Rotation[][3])
{
   double delta,gamma;                              
   //double Ry[3][3], Rz[3][3];

   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;      
   
   if (Vector[0] != 0)  
      delta = atan(Vector[2]/Vector[0]); 
   else 
      delta = pi/2.; 
   RotatesAboutY(delta,Vector); //,Ry); //## 01/11/27
   if (Vector[0] < 0)  
      gamma = 1/2.*atan(Vector[1]/Vector[0]);  
   else if (Vector[0] > 0)   {      
      if (atan(Vector[1]/Vector[0]) < 0) 
         gamma = 1/2.*(atan(Vector[1]/Vector[0]) + pi );
      else
         gamma = 1/2.*(atan(Vector[1]/Vector[0]) - pi );               
   }      
   else {
      if (Vector[1] > 0) 
         gamma = -pi/4;
      else
         gamma =  pi/4;       
   }      
   RotatesAboutZ(2*gamma,Vector); //,Rz); 

   thetaZthetaY[0] = 2*gamma; // ## Rz       
   thetaZthetaY[1] = delta;   // ## Ry   
        
   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;       
   //Multiply(Rotation,Rz,Ry);     
}


//## It produces the angles and the matrices that rotates the 
//#  vector to (0, y, z > 0)
//## By: Ivan Lima
void TransRyV3to_YZ( double *Vector, double *thetaY)
{
   double delta;                              

   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;      
   if (Vector[2] != 0)  {
      delta = -atan(Vector[0]/Vector[2]); 
   }   
   else {
      delta = -sgn(Vector[0])*pi/2.; 
   }   
   if (Vector[2] < 0 ) {
      delta += -pi;
    }  
      
   RotatesAboutY(delta,Vector); //,Ry);

   *thetaY = delta;   // ## Ry   
   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;       
   //Multiply(Rotation,Rz,Ry);     
}


//## Transforms from (vx,vy,vz) to (-1,0,0) 
//## Returning the angles and the Trans. Matrix RzRx
//## This method should be replaced by the one where the
//## two angles are separate scalar quantities (not a vector)
//## By: Ivan Lima
void TransRzRxV3to_Xhat( double *Vector, double *thetaZthetaX)
//                    ,double Rotation[][3])
{

   double thetaZ;
   double thetaX;
   TransRzRxV3to_Xhat( Vector, &thetaZ, &thetaX);
   thetaZthetaX[0] = thetaZ; // ## Rz angle       
   thetaZthetaX[1] = thetaX;   // ## Rx angle    
/*
   double delta,gamma;                                 
   
   //double Rx[3][3], Rz[3][3];

   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;      
   if (Vector[1] != 0)  
      delta = atan(Vector[2]/Vector[1]); 
   else 
      delta = pi/2.; 
   RotatesAboutX(delta,Vector); //,Rx);
   if (Vector[0] < 0)  
      gamma = 1/2.*atan(Vector[1]/Vector[0]);  
   else if (Vector[0] > 0)   {      
      if (atan(Vector[1]/Vector[0]) < 0) 
         gamma = 1/2.*(atan(Vector[1]/Vector[0]) + pi );
      else
         gamma = 1/2.*(atan(Vector[1]/Vector[0]) - pi );               
   }      
   else {
      if (Vector[1] > 0) 
         gamma = -pi/4;
      else
         gamma =  pi/4;       
   }      
   RotatesAboutZ(2*gamma,Vector); //,Rz); 

   thetaZthetaX[0] = 2*gamma; // ## Rz angle       
   thetaZthetaX[1] = delta;   // ## Rx angle    
   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;       
   //Multiply(Rotation,Rz,Rx);     
*/   
}   





//## Transforms from (vx,vy,vz) to (-1,0,0) 
//## Returning the angles and the Trans. Matrix RzRx
//## By: Ivan Lima
void TransRzRxV3to_Xhat( double *Vector, double *thetaZ, double *thetaX)
//                    ,double Rotation[][3])
{
   double delta,gamma;                              

   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;      
   if (Vector[1] != 0)  
      delta = atan(Vector[2]/Vector[1]); 
   else 
      delta = pi/2.; 
   RotatesAboutX(delta,Vector); //,Rx);
   if (Vector[0] < 0)  
      gamma = 1/2.*atan(Vector[1]/Vector[0]);  
   else if (Vector[0] > 0)   {      
      if (atan(Vector[1]/Vector[0]) < 0) 
         gamma = 1/2.*(atan(Vector[1]/Vector[0]) + pi );
      else
         gamma = 1/2.*(atan(Vector[1]/Vector[0]) - pi );               
   }      
   else {
      if (Vector[1] > 0) 
         gamma = -pi/4;
      else
         gamma =  pi/4;       
   }      
   RotatesAboutZ(2*gamma,Vector); //,Rz); 

   *thetaZ = 2*gamma; // ## Rz angle       
   *thetaX = delta;   // ## Rx angle    
   //cout<<"Trans "<<Vector[0]<<" "<<Vector[1]<<" "<<Vector[2]<<endl;       
   //Multiply(Rotation,Rz,Rx);     
}   









//## Solve the 3x3 multiplication: Right3x3 = Left3x3*Right3x3
//## By: Ivan Lima
void Multiply( double Left[][3], double Right[][3])
{
   double RightTmp[3][3];
   for (int ii = 0; ii < 3; ii++) {
      for (int jj = 0; jj < 3; jj++) {
         RightTmp[ii][jj] = Right[ii][jj];             
      }   
   }   
   for (int ii = 0; ii < 3; ii++) {
      for (int jj = 0; jj < 3; jj++) {
         for (int kk = 0; kk < 3; kk++) {
            Right[ii][jj] += Left[ii][kk]*RightTmp[kk][jj];             
         }   
         //cout << "  "<<  ii << " " << jj << "  " << Left[ii][jj] << endl;          
         //cout << "  "<<  ii << " " << jj << "  " << Right[ii][jj] << endl;          
         //cout << ii << " " << jj << "  " << Result[ii][jj] << endl;          
      }   
   }   
}


//## Solve the matrix3x3 vector multipl.: Right3 = Left3x3*Right3
//## By: Ivan Lima
void Multiply(double Left[][3],double Right[3])
{
   double RightTmp[3];
   for (int ii = 0; ii < 3; ii++)
      RightTmp[ii] = Right[ii];
   for (int ii = 0; ii < 3; ii++) {   
      Right[ii] = 0;
      for (int jj = 0; jj < 3; jj++) {
         Right[ii] += Left[ii][jj]*RightTmp[jj];             
      }   
  }       
}

//## Make the CrossProduct

//## The result will be returned in Right
//## By: Ivan Lima
void CrossProduct(double *Left, double *Right, int FirstCoord)
{
   //## FirstCoord determines where the vector starts
   //#  0 -> 0,1,2
   //#  1 -> 1,2,3
   double RightTmp[4];
   for (int ii = 0; ii < 3; ii++) {
      RightTmp[ii+FirstCoord] = Right[ii+FirstCoord];
   }
   Right[0+FirstCoord] = Left[1+FirstCoord]*RightTmp[2+FirstCoord] 
                      -Left[2+FirstCoord]*RightTmp[1+FirstCoord];
   Right[1+FirstCoord] = Left[2+FirstCoord]*RightTmp[0+FirstCoord] 
                      -Left[0+FirstCoord]*RightTmp[2+FirstCoord];
   Right[2+FirstCoord] = Left[0+FirstCoord]*RightTmp[1+FirstCoord] 
                      -Left[1+FirstCoord]*RightTmp[0+FirstCoord];   
}


//## Rotates about Z
//## By: Ivan Lima
void RotatesAboutZ(double angle, double *Vector) //, double Rz[][3])
{

   double Rz[3][3];
   //for (int ii=0; ii < 3; ii++)
   //   for (int jj = 0; jj < 3; jj++)
   //      Rz[ii][jj] = 0;    
   Rz[0][0] =  cos(angle);   Rz[0][1] =  sin(angle);
   Rz[1][0] = -sin(angle);   Rz[1][1] =  cos(angle);
   //Rz[2][2] = 1.;
   
   double tmpVector[3];
   tmpVector[0] = Vector[0];
   tmpVector[1] = Vector[1];    
   Vector[0] = Rz[0][0]*tmpVector[0] + Rz[0][1]*tmpVector[1];
   Vector[1] = Rz[1][0]*tmpVector[0] + Rz[1][1]*tmpVector[1];      
}


//## Rotates about X
//## By: Ivan Lima
void RotatesAboutX(double angle, double *Vector) //, double Rx[][3])
{    
   double Rx[3][3];
   //for (int ii = 0; ii < 3; ii++)
   //   for (int jj=0; jj < 3; jj++)
   //      Rx[ii][jj] = 0;    
   //Rx[0][0] = 1.;         
   Rx[1][1] =  cos(angle);    Rx[1][2] =  sin(angle);
   Rx[2][1] = -sin(angle);    Rx[2][2] =  cos(angle);
   
   double tmpVector[3];
   tmpVector[1] = Vector[1];   
   tmpVector[2] = Vector[2];   
   Vector[1] = Rx[1][1]*tmpVector[1] + Rx[1][2]*tmpVector[2];
   Vector[2] = Rx[2][1]*tmpVector[1] + Rx[2][2]*tmpVector[2];   
}


// ## Rotates about Y
//##  By: Ivan Lima
void RotatesAboutY(double angle, double *Vector) //, double Ry[][3])
{
   double Ry[3][3];
   //for (int ii = 0; ii < 3; ii++)
   //   for (int jj=0; jj < 3; jj++)
   //      Ry[ii][jj] = 0;  
   Ry[0][0] =  cos(angle);              Ry[0][2] =  sin(angle);           
   //                        Ry[1][1] = 1.;         
   Ry[2][0] = -sin(angle);               Ry[2][2] =  cos(angle);
   
   double tmpVector[3];
   tmpVector[0] = Vector[0];   
   tmpVector[2] = Vector[2];   
   Vector[0] = Ry[0][0]*tmpVector[0] + Ry[0][2]*tmpVector[2];
   Vector[2] = Ry[2][0]*tmpVector[0] + Ry[2][2]*tmpVector[2];   
}


// ## Scalar product
//## The coordinates are from 0 to 2
//## By: Ivan Lima
double ScalarProduct(double *VectorA, double *VectorB) // from 0 to 2
{
   return( VectorA[0]*VectorB[0]
          +VectorA[1]*VectorB[1]
          +VectorA[2]*VectorB[2]);
} 

//## Gets the length of a vector
//## By: Ivan Lima
double GetVectorLength(double *Vector, int Dimension) // from 0 to dimension-1
{
   double VectorLength = 0;
   for (int ii = 0; ii < Dimension; ii++) {
      VectorLength += sq(Vector[ii]);
   }   
   return sqrt(VectorLength);   
} 



//## Rotates a Jones vector about X
//## By: Ivan Lima
void RotatesAboutX(double angle, cplx *JonesVector) 
{    
   cplx TransM[2][2];     
   TransM[0][0] =  exp(-jc*angle/2.);  TransM[0][1] =  0;
   TransM[1][0] = 0;                  TransM[1][1] =  exp(jc*angle/2.);
   
   cplx tmpJonesVector[2];
   tmpJonesVector[0] = JonesVector[0];   
   tmpJonesVector[1] = JonesVector[1];   
   JonesVector[0] =  TransM[0][0]*tmpJonesVector[0] 
                   + TransM[0][1]*tmpJonesVector[1];
   JonesVector[1] =  TransM[1][0]*tmpJonesVector[0] 
                   + TransM[1][1]*tmpJonesVector[1];   
}


//## Rotates a Jones vector about Y
//## By: Ivan Lima
void RotatesAboutY(double angle, cplx *JonesVector) 
{    
   cplx TransM[2][2];     
   TransM[0][0] =  cos(angle/2.);    TransM[0][1] =  jc*sin(angle/2.);
   TransM[1][0] = jc*sin(angle/2.); TransM[1][1] =  cos(angle/2.);
   
   cplx tmpJonesVector[2];
   tmpJonesVector[0] = JonesVector[0];   
   tmpJonesVector[1] = JonesVector[1];   
   JonesVector[0] =  TransM[0][0]*tmpJonesVector[0] 
                   + TransM[0][1]*tmpJonesVector[1];
   JonesVector[1] =  TransM[1][0]*tmpJonesVector[0] 
                   + TransM[1][1]*tmpJonesVector[1];   
}




//## Rotates a Jones vector about Z
//## By: Ivan Lima
void RotatesAboutZ(double angle, cplx *JonesVector) 
{    
   cplx TransM[2][2];     
   TransM[0][0] =  cos(angle/2.);    TransM[0][1] =  sin(angle/2.);
   TransM[1][0] = -sin(angle/2.);    TransM[1][1] =  cos(angle/2.);
   
   cplx tmpJonesVector[2];
   tmpJonesVector[0] = JonesVector[0];   
   tmpJonesVector[1] = JonesVector[1];   
   JonesVector[0] =  TransM[0][0]*tmpJonesVector[0] 
                   + TransM[0][1]*tmpJonesVector[1];
   JonesVector[1] =  TransM[1][0]*tmpJonesVector[0] 
                   + TransM[1][1]*tmpJonesVector[1];   
}




// ############################################
 
//## Converts from dB to the linear scale
//## By: Ivan Lima
double dB2Linear(double dBvalue)
{
   return(pow(10.,dBvalue/10.) ); 
}


// ########################################

double dBm2Linear(double dBmPower)
{

  // Written by John Zweck, Dec 15th 2000

  // Returns power is in Watts, W

  return(pow(10.0,dBmPower/10.0 - 3.0) );


}

// ####################################

double Linear2dB(double LinearValue)
{
 // Written by John Zweck, Dec 15th 2000

  return 10.0*log10(LinearValue);  


}

// #####################################

double Linear2dBm(double LinearPower)
{
 // Written by John Zweck, Dec 15th 2000

  // The input power LinearValue must be power in Watts, W.

   return 10.0*log10(LinearPower*1000.0);  

}


//####################################
//## This function is adequate to limit
//#  the maximum amount of power penalty
//#  that a system can suffer
//## By: Ivan Lima (5/10/02)
double Linear2dB_Limited(double LinearValue, double MaxValue_dB)
{
  //if (LinearValue <= 0) {
  if ( LinearValue >= dB2Linear(MaxValue_dB) ) {
     //## Maximum value set to 20 dB
     LinearValue = dB2Linear(MaxValue_dB);
  }   
  return Linear2dB(LinearValue);   
}


//#####################################
//## By: Ivan Lima
double GetDegreeOfPolarization(double *StokesVector)
{
   if ( StokesVector[0]*(1 + 2*NEAR_ZERO) < 
         sqrt(sq(StokesVector[1]) + sq(StokesVector[2])
        +sq(StokesVector[3])) ) {
      cout << "From: DegreeOfPolarization()" << endl 
           << "Error: The  Stokes  Vector has DOP > 1." 
           << endl
           << "Please, check simulation algorithm" 
           << endl;
      for (int ii = 0; ii < 4; ii++)
         cout << "StokesVector[" << ii << "] = " 
              << StokesVector[ii] << endl;     
      exit(1);
   }   
   return(sqrt( sq(StokesVector[1]) + sq(StokesVector[2])
               +sq(StokesVector[3]) )/StokesVector[0] );
}      
     
//######################################################     
//## By: Ivan Lima
double GetPhase(const cplx Phasor) 
//## Garanteed from -pi to pi
{
   if (abs(Phasor) == 0) 
      return(0.);
   else   // phase = sgn( Aim ) * acos( Are/Ao)
      return( sgn(real(-jc*Phasor))*acos(real(Phasor)/abs(Phasor)) );      
}


//##################################################
//## By: Ivan Lima
void Stokes2Jones(const double *Stokes, cplx *Jones)
//## The Stokes vector should have dimension 4: from 1 to 3   
//## The Jones  vector should have dimension 2: from 0 to 1  
//## The Jones  vector is returned normalized
{   
   //## Ths step provides robustness agains lack of normalization
   //#  Obviously, it assumes that the Stokes parameters
   //#  are totaly polarized.
   double Stokes0 = sqrt(sq(Stokes[1])+sq(Stokes[2])+sq(Stokes[3]) );
   //if (Stokes0 < NEAR_ZERO) {
   //   cout << "Error: the Stokes parameters cannot be zero." << endl;
   //   cout << "Exiting from Stokes2Jones()" << endl;
   //   exit(1);
   //}             
   if (Stokes0 == 0) {
      Jones[0] = Jones[1] = 0;
   }
   else {
      double StokesNorm[4];   
      StokesNorm[1] = Stokes[1]/Stokes0;
      StokesNorm[2] = Stokes[2]/Stokes0;
      StokesNorm[3] = Stokes[3]/Stokes0;   
      double PhaseXY; 
      //if (sq(StokesNorm[1]) >= 1 - NEAR_ZERO) {
      if (sq(StokesNorm[1]) == 1) {      
            PhaseXY = 0;
      }      
      else {
         PhaseXY = acos( StokesNorm[2]/sqrt(1-sq(StokesNorm[1]) ) );
      }   
      Jones[0] = sqrt((1. + StokesNorm[1])/2.);
      Jones[1] = sqrt((1. - StokesNorm[1])/2.)
                *exp(-jc*sgn(StokesNorm[3])*PhaseXY);         
   }          
}

//#####################################################
//## By: Ivan Lima
void Stokes2JonesKeepPower(const double *Stokes, cplx *Jones)
{
   Stokes2Jones(Stokes,Jones);
   double TotPower = sqrt(sqrt( sq(Stokes[1])
                          +sq(Stokes[2])+sq(Stokes[3])) );
   Jones[0] *= TotPower; 
   Jones[1] *= TotPower;                      
}


//##################################################
//## By: Ivan Lima
void Jones2Stokes(const cplx *Jones, double *Stokes)
//## The Jones  vector should have dimension 2: from 0 to 1  
//## The Stokes vector should have dimension 4: from 0 to 3 
//## The Stokes vector is returned normalized
{   
   Stokes[0]   = sq(abs(Jones[0])) + sq(abs(Jones[1]));
   if (Stokes[0] != 0) {
      Stokes[1] = (sq(abs(Jones[0]))-sq(abs(Jones[1])) )
                  /Stokes[0];                           //##  |uc|^2-|vc|^2
      Stokes[2] = 2*real(Jones[0]*conj(Jones[1]))
                  /Stokes[0];                           //##  2 Re{uc . vc*}
      Stokes[3] = 2*real(-jc*Jones[0]*conj(Jones[1]))  
                  /Stokes[0];                           //## 2 Im{uc . vc*}         
      Stokes[0] = 1;             
   }
   else  {
      Stokes[1] = Stokes[2] = Stokes[3] = 0;
   }                     
}


//##################################################
//## By: Ivan Lima
void Jones2StokesKeepPower(const cplx *Jones, double *Stokes)
{
   Jones2Stokes(Jones,Stokes);
   double TotPower = sq(abs(Jones[0])) + sq(abs(Jones[1]));
   for (int ii = 0; ii < 4; ii++) {
      Stokes[ii] *= TotPower;
   }
}



// ##### End of functions written by Ivan T. Lima Jr.


// #################################################################
// #################################################################

// Functions written by John Zweck

double Frequency2Wavelength(double Frequency)
{
  // Frequency is in Hz

  return LightSpeed/Frequency;  
}

// ############################################################

double Wavelength2Frequency(double Wavelength)
{
  return LightSpeed/Wavelength; // returns Frequency  in Hz
}


// ##############################################################
// ##############################################################

// ##### The functions below were written by Ronald Holzloehner

//--------------------------------------------
// Box-Mueller random number generator
// see D. Knuth: "The Art of Computer Programming", 3rd ed, vol 2, p.122

void gaussian_pdf(int N, double stan_dev, double *vec)
{
  double *p_vec = &(vec[0]);
  double *const p_vec_end = &(vec[N]);
  double v1, v2, SSs;

  for( ;p_vec < p_vec_end; )
  {
    do{ 
        v1 = 2.0*drand48() - 1.0;
        v2 = 2.0*drand48() - 1.0;
        SSs = sq(v1) + sq(v2);
      }while(SSs > 1.0);

    double a = stan_dev * sqrt(-2.0 * log(SSs)/SSs);
    *p_vec = a * v1;
    p_vec++;
    *p_vec = a * v2;
    p_vec++;
  }
}


//--------------------------------------------------------------------
// 5th order normalized Bessel function, see bessel5.m

cplx bessel5(double omega, double omega_0, double FWHM)
{
  // 5th order
  const double mult = 4.8548214044214655;

  const cplx s = cplx( 0.0, (omega-omega_0)/FWHM * mult );

  // 5th order normalized polynomial
  return 945.0  /  (945.0 + s * (945.0 + s * (420.0+s * (105.0+s * (15.0+s)  ))));
}


//----------------------------------------------------------------------
// electrical filter, 5th-order Bessel
// Input: my_power_time[]
// Output: my_power_time[]
// global side effects: none

void electrical_filter(cplx *my_power_time, cfftw *fft, double fwhm_norm)
{
  fft->FFT(my_power_time, 0);
  cplx *power_f_p = fft->work;
  const int N = fft->give_N();
  double phase0125 = arg(bessel5(fwhm_norm/8.0, 0.0, fwhm_norm)); // go almost 0.5*fwhm
  for(int i=0; i < N; i++, power_f_p++)
  {
    // delta_f is the same for the splitstep data vector (length N_SS)
    // and for downsampled values b[] since delta_f = 1/T_total
    // CAUTION: delta_f is independent of step number, delta_t not!!
    double freq = fft->freq(i);
    double phase = -freq/(fwhm_norm/8.0) * phase0125; // minus to compensate
    *power_f_p *= cplx(cos(phase), sin(phase))
      * bessel5(freq, 0.0, fwhm_norm);
  }
  // transform back to time domain
  fft->IFFT(0, my_power_time);
  
//    cout << "\n  Values at cutoff (should be 1/sqrt(2)): " << 
//    abs(bessel5(bessel_FWHM/2, 0, bessel_FWHM)) << "\t" <<
//    abs(bessel5(-bessel_FWHM/2, 0, bessel_FWHM)) << endl;
} // electrical_filter()


//----------------------------------------------------------------------
// implements optical Gaussian filter
// multiplies optical field (not intensity!) in frequency domain
// Input: my_field_freq[]
// Output: both  my_field_time[], my_field_freq[]
// CAUTION: gives FWHM in intensity, not field !! (checked 21 Jan 00)

void optical_filter(cplx *my_field_time, cplx *my_field_freq, 
   cfftw *fft, double fwhm_norm)
{
  cplx *field_f = my_field_freq;
  int i;
  double sigma = fwhm_norm / 2.35482 * M_SQRT2;
  const int N = fft->give_N();
  for(i=0; i < N; i++, field_f++)
  {
    // CAUTION: delta_f is independent of step number, delta_t not!!
    *field_f *= exp(-0.5*sq(fft->freq(i)/sigma));
  }
  fft->IFFT(my_field_freq, my_field_time);

} // optical_filter()

//-----------------------------------------------------------------------
// enable program to be migrated by MOSIX
// Options: MOSIX_UNLOCK or MOSIX_LOCK
int lock_or_unlock_mosix(int what)
{
  ofstream *lock_file;

  char lock_name[50];
  sprintf(lock_name, "/proc/%d/lock", getpid());
  lock_file = new ofstream( lock_name, ios::out );
  if(lock_file->bad())
    return 1;
  *lock_file << ((what==MOSIX_LOCK)?"1":"0") << endl;
  lock_file->close();
  return 0;
}


//-----------------------------------------------------------------------
// tell MOSIX this is a pure CPU job (as opposed to one with lots of I/O)

int cpujob_mosix(void)
{
  ofstream *cpu_file;

  cpu_file = new ofstream( "/proc/mosix/decay/cpujob", ios::out );
  if(cpu_file->bad())
    return 1;
  *cpu_file << 1 << endl;
  cpu_file->close();
  return 0;
}


//-----------------------------------------------------------------------
// tell MOSIX that slow statistical decay applies to this program
int slow_mosix(void)
{
  ofstream *slow_file;

  slow_file = new ofstream( "/proc/mosix/decay/slow", ios::out );
  if(slow_file->bad())
    return 1;
  *slow_file << 1 << endl;
  slow_file->close();
  return 0;
}


//-----------------------------------------------------------------------
double peak_power(cplx *in, int N)
{
  cplx *p = in;
  cplx *const p_end = in+N;
  double max_power = norm(*p);

  for(p++; p < p_end; p++)
  {
    double tmp = norm(*p);
    if(tmp > max_power)
      max_power = tmp;
  }
  return max_power;
}


//-----------------------------------------------------------------------
// Simpson's rule with a little trick -- it is assumed that
// the signal is periodic and in[0] == in[N]
// See Springer's Mathe, p. 377
//
// *in: input array
// N:   length of array
// width: 'time window', real step width times N

double pulse_energy(cplx *in, int N, double width)
{
  if(N%2)
  {
    cout << "\a\npulse_energy(): N must be even! Returning NAN..." << endl;
    cout << "Error in pulse_energy()." << endl;
    exit(1);
    //return NAN;
  }

  double sum1 = 0.0;
  cplx *p = in+1;
  cplx *const p_end = in+N;
  for(; p < p_end; p+=2)
    sum1 += norm(*p);

  double sum2 = 0.0;
  p = in; // yes, we start from index 0 here -- this inludes in[0]
  for(; p < p_end; p+=2)
    sum2 += norm(*p);

  // values in[0] and in[N] = in[0] are included in sum2
  return  width/double(N)/3.0 * ( 4.0*sum1 + 2.0*sum2 );
}
// ##### End functions written by Ronald Holzloehner


// ##### Functions from Numerical Recipes

// ## Attention: all the random numbers inside the program are 
// obtained calling ranG(...). ranG(...) is defined here as
// ran2(..), from Numerical Recipes, by default. 
// But it can be changed to any other function alike.

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

//float ran2(long *idum)  //##### Original name from Numerical Recipes
float ranG(long *idum)
{
        int j;
        long k;
        static long idum2=123456789;
        static long iy=0;
        static long iv[NTAB];
        float temp;

        if (*idum <= 0) {
                if (-(*idum) < 1) *idum=1;
                else *idum = -(*idum);
                idum2=(*idum);
                for (j=NTAB+7;j>=0;j--) {
                        k=(*idum)/IQ1;
                        *idum=IA1*(*idum-k*IQ1)-k*IR1;
                        if (*idum < 0) *idum += IM1;
                        if (j < NTAB) iv[j] = *idum;
                }
                iy=iv[0];
        }
        k=(*idum)/IQ1;
        *idum=IA1*(*idum-k*IQ1)-k*IR1;
        if (*idum < 0) *idum += IM1;
        k=idum2/IQ2;
        idum2=IA2*(idum2-k*IQ2)-k*IR2;
        if (idum2 < 0) idum2 += IM2;
        j=iy/NDIV;
        iy=iv[j]-idum2;
        iv[j] = *idum;
        if (iy < 1) iy += IMM1;        
        if ((temp=AM*iy) > RNMX) return RNMX;
        else return temp;
}
#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!#:V3+). 


#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran1(long *idum)  //##### Original name from Numerical Recipes
//float ranG(long *idum)
{
        int j;
        long k;
        static long iy=0;
        static long iv[NTAB];
        float temp;

        if (*idum <= 0 || !iy) {
                if (-(*idum) < 1) *idum=1;
                else *idum = -(*idum);
                for (j=NTAB+7;j>=0;j--) {
                        k=(*idum)/IQ;
                        *idum=IA*(*idum-k*IQ)-IR*k;
                        if (*idum < 0) *idum += IM;
                        if (j < NTAB) iv[j] = *idum;
                }
                iy=iv[0];
        }
        k=(*idum)/IQ;
        *idum=IA*(*idum-k*IQ)-IR*k;
        if (*idum < 0) *idum += IM;
        j=iy/NDIV;
        iy=iv[j];
        iv[j] = *idum;
        if ((temp=AM*iy) > RNMX) return RNMX;
        else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!#:V3+). 


// ## Changed by Ivan in 8/25/2000 respect to the random generator
// #  For arbitrary mean and standart deviation use it like:
// #  myRandNumber = myMean + myStdDev*Gaussian_pdf(mySeedPointer);
#include <math.h>

float Gaussian_pdf(long *idum)
//float gasdev(long *idum)
{
        //float ran1(long *idum);
        float ranG(long *idum);
        static int iset=0;
        static float gset;
        float fac,rsq,v1,v2;

        if  (iset == 0) {
                do {
                        v1=2.0*ranG(idum)-1.0;
                        v2=2.0*ranG(idum)-1.0;
                        //v1=2.0*ran1(idum)-1.0;
                        //v2=2.0*ran1(idum)-1.0;                        
                        rsq=v1*v1+v2*v2;
                } while (rsq >= 1.0 || rsq == 0.0);
                fac=sqrt(-2.0*log(rsq)/rsq);
                gset=v1*fac;
                iset=1;
                return v2*fac;
        } else {
                iset=0;
                return gset;
        }
}
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+).


// ## Optimization functions from Numerical Recipes
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 200

void powell(float p[], float xi[][4], int n, float ftol, int *iter, float *fret,
        float (*func)(float []))
//void powell(float p[], float **xi, int n, float ftol, int *iter, float *fret,
//      float (*func)(float []))        
{
        void linmin(float p[], float xi[], int n, float *fret,
                float (*func)(float []));
        int i,ibig,j;
        float del,fp,fptt,t,*pt,*ptt,*xit;

        pt=nrvector(1,n);
        ptt=nrvector(1,n);
        xit=nrvector(1,n);
        *fret=(*func)(p);
        for (j=1;j<=n;j++) pt[j]=p[j];
        for (*iter=1;;++(*iter)) {
                fp=(*fret);
                ibig=0;
                del=0.0;
                for (i=1;i<=n;i++) {
                        for (j=1;j<=n;j++) xit[j]=xi[j][i];
                        fptt=(*fret);
                        linmin(p,xit,n,fret,func);
                        if (fabs(fptt-(*fret)) > del) {
                                del=fabs(fptt-(*fret));
                                ibig=i;
                        }
                }
                if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
                        free_nrvector(xit,1,n);
                        free_nrvector(ptt,1,n);
                        free_nrvector(pt,1,n);
                        return;
                }
                if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");
                for (j=1;j<=n;j++) {
                        ptt[j]=2.0*p[j]-pt[j];
                        xit[j]=p[j]-pt[j];
                        pt[j]=p[j];
                }
                fptt=(*func)(ptt);
                if (fptt < fp) {
                        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
                        if (t < 0.0) {
                                linmin(p,xit,n,fret,func);
                                for (j=1;j<=n;j++) {
                                        xi[j][ibig]=xi[j][n];
                                        xi[j][n]=xit[j];
                                }
                        }
                }
        }
}
#undef ITMAX
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 


#define NRANSI
#include "nrutil.h"
#define TOL 2.0e-4

int ncom;
float *pcom,*xicom,(*nrfunc)(float []);

void linmin(float p[], float xi[], int n, float *fret, float (*func)(float []))
{
        float brent(float ax, float bx, float cx,
                float (*f)(float), float tol, float *xmin);
        float f1dim(float x);
        void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,
                float *fc, float (*func)(float));
        int j;
        float xx,xmin,fx,fb,fa,bx,ax;

        ncom=n;
        pcom=nrvector(1,n);
        xicom=nrvector(1,n);
        nrfunc=func;
        for (j=1;j<=n;j++) {
                pcom[j]=p[j];
                xicom[j]=xi[j];
        }
        ax=0.0;
        xx=1.0;
        mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
        *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);        
        for (j=1;j<=n;j++) {
                xi[j] *= xmin;
                p[j] += xi[j];
        }
        free_nrvector(xicom,1,n);
        free_nrvector(pcom,1,n);
}
#undef TOL
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 

#define NRANSI
#include "nrutil.h"

extern int ncom;
extern float *pcom,*xicom,(*nrfunc)(float []);

float f1dim(float x)
{
        int j;
        float f,*xt;

        xt=nrvector(1,ncom);
        for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
        f=(*nrfunc)(xt);
        free_nrvector(xt,1,ncom);
        return f;
}
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 

#include <math.h>
#define NRANSI
#include "nrutil.h"
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY 1.0e-20
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb, float *fc,
        float (*func)(float))
{
        float ulim,u,r,q,fu,dum;

        *fa=(*func)(*ax);
        *fb=(*func)(*bx);
        if (*fb > *fa) {
                SHFT(dum,*ax,*bx,dum)
                SHFT(dum,*fb,*fa,dum)
        }
        *cx=(*bx)+GOLD*(*bx-*ax);
        *fc=(*func)(*cx);
        while (*fb > *fc) {
                r=(*bx-*ax)*(*fb-*fc);
                q=(*bx-*cx)*(*fb-*fa);
                u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
                        (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
                ulim=(*bx)+GLIMIT*(*cx-*bx);
                if ((*bx-u)*(u-*cx) > 0.0) {
                        fu=(*func)(u);
                        if (fu < *fc) {
                                *ax=(*bx);
                                *bx=u;
                                *fa=(*fb);
                                *fb=fu;
                                return;
                        } else if (fu > *fb) {
                                *cx=u;
                                *fc=fu;
                                return;
                        }
                        u=(*cx)+GOLD*(*cx-*bx);
                        fu=(*func)(u);
                } else if ((*cx-u)*(u-ulim) > 0.0) {
                        fu=(*func)(u);
                        if (fu < *fc) {
                                SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
                                SHFT(*fb,*fc,fu,(*func)(u))
                        }
                } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
                        u=ulim;
                        fu=(*func)(u);
                } else {
                        u=(*cx)+GOLD*(*cx-*bx);
                        fu=(*func)(u);
                }
                SHFT(*ax,*bx,*cx,u)
                SHFT(*fa,*fb,*fc,fu)
        }
}
#undef GOLD
#undef GLIMIT
#undef TINY
#undef SHFT
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 

#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 100
#define CGOLD 0.3819660
#define ZEPS 1.0e-10
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

float brent(float ax, float bx, float cx, float (*f)(float), float tol,
        float *xmin)
{
        int iter;
        float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
        float e=0.0;

        a=(ax < cx ? ax : cx);
        b=(ax > cx ? ax : cx);
        x=w=v=bx;
        fw=fv=fx=(*f)(x);
        for (iter=1;iter<=ITMAX;iter++) {
                xm=0.5*(a+b);
                tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
                if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
                        *xmin=x;
                        return fx;
                }
                if (fabs(e) > tol1) {
                        r=(x-w)*(fx-fv);
                        q=(x-v)*(fx-fw);
                        p=(x-v)*q-(x-w)*r;
                        q=2.0*(q-r);
                        if (q > 0.0) p = -p;
                        q=fabs(q);
                        etemp=e;
                        e=d;
                        if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
                                d=CGOLD*(e=(x >= xm ? a-x : b-x));
                        else {
                                d=p/q;
                                u=x+d;
                                if (u-a < tol2 || b-u < tol2)
                                        d=SIGN(tol1,xm-x);
                        }
                } else {
                        d=CGOLD*(e=(x >= xm ? a-x : b-x));
                }
                u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
                fu=(*f)(u);
                if (fu <= fx) {
                        if (u >= x) a=x; else b=x;
                        SHFT(v,w,x,u)
                        SHFT(fv,fw,fx,fu)
                } else {
                        if (u < x) a=u; else b=u;
                        if (fu <= fw || w == x) {
                                v=w;
                                w=u;
                                fv=fw;
                                fw=fu;
                        } else if (fu <= fv || v == x || v == w) {
                                v=u;
                                fv=fu;
                        }
                }
        }
        nrerror("Too many iterations in brent");
        *xmin=x;
        return fx;
}
#undef ITMAX
#undef CGOLD
#undef ZEPS
#undef SHFT
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 




// ## Conjugate Gradient Method in double precision
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 200
#define EPS 1.0e-10
#define FREEALLD free_dvector(xi,1,n);free_dvector(h,1,n);free_dvector(g,1,n);

void frprmnd(double p[], int n, double ftol, int *iter, double *fret,
        double (*func)(double []), void (*dfunc)(double [], double []))
{
        void linmind(double p[], double xi[], int n, double *fret,
                double (*func)(double []));
        int j,its;
        double gg,gam,fp,dgg;
        double *g,*h,*xi;

        g=dvector(1,n);
        h=dvector(1,n);
        xi=dvector(1,n);
        fp=(*func)(p);
        (*dfunc)(p,xi);
        for (j=1;j<=n;j++) {
                g[j] = -xi[j];
                xi[j]=h[j]=g[j];
        }
        for (its=1;its<=ITMAX;its++) {
                *iter=its;
                linmind(p,xi,n,fret,func);
                if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {
                        FREEALLD
                        return;
                }
                fp=(*func)(p);
                (*dfunc)(p,xi);
                dgg=gg=0.0;
                for (j=1;j<=n;j++) {
                        gg += g[j]*g[j];
                        dgg += (xi[j]+g[j])*xi[j];
                }
                if (gg == 0.0) {
                        FREEALLD
                        return;
                }
                gam=dgg/gg;
                for (j=1;j<=n;j++) {
                        g[j] = -xi[j];
                        xi[j]=h[j]=g[j]+gam*h[j];
                }
        }
        nrerror("Too many iterations in frprmn");
}
#undef ITMAX
#undef EPS
#undef FREEALLD
#undef NRANSI
// (C) Copr. 19 


#define NRANSI
#include "nrutil.h"
#define TOL 2.0e-4

//int ncom;
double *pcomd,*xicomd,(*nrfuncd)(double []);

void linmind(double p[], double xi[], int n, double *fret, double (*func)(double []))
{
        double brentd(double ax, double bx, double cx,
                double (*f)(double), double tol, double *xmin);
        double f1dim(double x);
        void mnbrakd(double *ax, double *bx, double *cx, double *fa, double *fb,
                double *fc, double (*func)(double));
        int j;
        double xx,xmin,fx,fb,fa,bx,ax;

        ncom=n;
        pcomd=dvector(1,n);
        xicomd=dvector(1,n);
        nrfuncd=func;
        for (j=1;j<=n;j++) {
                pcomd[j]=p[j];
                xicomd[j]=xi[j];
        }
        ax=0.0;
        xx=1.0;
        mnbrakd(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
        *fret=brentd(ax,xx,bx,f1dim,TOL,&xmin);        
        for (j=1;j<=n;j++) {
                xi[j] *= xmin;
                p[j] += xi[j];
        }
        free_dvector(xicomd,1,n);
        free_dvector(pcomd,1,n);
}
#undef TOL
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 

#define NRANSI
#include "nrutil.h"

extern int ncom;
extern double *pcomd,*xicomd,(*nrfuncd)(double []);

double f1dim(double x)
{
        int j;
        double f,*xt;

        xt=dvector(1,ncom);
        for (j=1;j<=ncom;j++) xt[j]=pcomd[j]+x*xicomd[j];
        f=(*nrfuncd)(xt);
        free_dvector(xt,1,ncom);
        return f;
}
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 

#include <math.h>
#define NRANSI
#include "nrutil.h"
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY 1.0e-20
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

void mnbrakd(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
        double (*func)(double))
{
        double ulim,u,r,q,fu,dum;

        *fa=(*func)(*ax);
        *fb=(*func)(*bx);
        if (*fb > *fa) {
                SHFT(dum,*ax,*bx,dum)
                SHFT(dum,*fb,*fa,dum)
        }
        *cx=(*bx)+GOLD*(*bx-*ax);
        *fc=(*func)(*cx);
        while (*fb > *fc) {
                r=(*bx-*ax)*(*fb-*fc);
                q=(*bx-*cx)*(*fb-*fa);
                u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
                        (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
                ulim=(*bx)+GLIMIT*(*cx-*bx);
                if ((*bx-u)*(u-*cx) > 0.0) {
                        fu=(*func)(u);
                        if (fu < *fc) {
                                *ax=(*bx);
                                *bx=u;
                                *fa=(*fb);
                                *fb=fu;
                                return;
                        } else if (fu > *fb) {
                                *cx=u;
                                *fc=fu;
                                return;
                        }
                        u=(*cx)+GOLD*(*cx-*bx);
                        fu=(*func)(u);
                } else if ((*cx-u)*(u-ulim) > 0.0) {
                        fu=(*func)(u);
                        if (fu < *fc) {
                                SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
                                SHFT(*fb,*fc,fu,(*func)(u))
                        }
                } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
                        u=ulim;
                        fu=(*func)(u);
                } else {
                        u=(*cx)+GOLD*(*cx-*bx);
                        fu=(*func)(u);
                }
                SHFT(*ax,*bx,*cx,u)
                SHFT(*fa,*fb,*fc,fu)
        }
}
#undef GOLD
#undef GLIMIT
#undef TINY
#undef SHFT
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 

#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 100
#define CGOLD 0.3819660
#define ZEPS 1.0e-10
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

double brentd(double ax, double bx, double cx, double (*f)(double), double tol,
        double *xmin)
{
        int iter;
        double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
        double e=0.0;

        a=(ax < cx ? ax : cx);
        b=(ax > cx ? ax : cx);
        x=w=v=bx;
        fw=fv=fx=(*f)(x);
        for (iter=1;iter<=ITMAX;iter++) {
                xm=0.5*(a+b);
                tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
                if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
                        *xmin=x;
                        return fx;
                }
                if (fabs(e) > tol1) {
                        r=(x-w)*(fx-fv);
                        q=(x-v)*(fx-fw);
                        p=(x-v)*q-(x-w)*r;
                        q=2.0*(q-r);
                        if (q > 0.0) p = -p;
                        q=fabs(q);
                        etemp=e;
                        e=d;
                        if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
                                d=CGOLD*(e=(x >= xm ? a-x : b-x));
                        else {
                                d=p/q;
                                u=x+d;
                                if (u-a < tol2 || b-u < tol2)
                                        d=SIGN(tol1,xm-x);
                        }
                } else {
                        d=CGOLD*(e=(x >= xm ? a-x : b-x));
                }
                u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
                fu=(*f)(u);
                if (fu <= fx) {
                        if (u >= x) a=x; else b=x;
                        SHFT(v,w,x,u)
                        SHFT(fv,fw,fx,fu)
                } else {
                        if (u < x) a=u; else b=u;
                        if (fu <= fw || w == x) {
                                v=w;
                                w=u;
                                fv=fw;
                                fw=fu;
                        } else if (fu <= fv || v == x || v == w) {
                                v=u;
                                fv=fu;
                        }
                }
        }
        nrerror("Too many iterations in brentd");
        *xmin=x;
        return fx;
}
#undef ITMAX
#undef CGOLD
#undef ZEPS
#undef SHFT
#undef NRANSI
// (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). 


// ## End of functions for the optimization in double precision




// ## Conjugate Gradient Method
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 200
#define EPS 1.0e-10
#define FREEALL free_nrvector(xi,1,n);free_nrvector(h,1,n);free_nrvector(g,1,n);

void frprmn(float p[], int n, float ftol, int *iter, float *fret,
        float (*func)(float []), void (*dfunc)(float [], float []))
{
        void linmin(float p[], float xi[], int n, float *fret,
                float (*func)(float []));
        int j,its;
        float gg,gam,fp,dgg;
        float *g,*h,*xi;

        g=nrvector(1,n);
        h=nrvector(1,n);
        xi=nrvector(1,n);
        fp=(*func)(p);
        (*dfunc)(p,xi);
        for (j=1;j<=n;j++) {
                g[j] = -xi[j];
                xi[j]=h[j]=g[j];
        }
        for (its=1;its<=ITMAX;its++) {
                *iter=its;
                linmin(p,xi,n,fret,func);
                if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {
                        FREEALL
                        return;
                }
                fp=(*func)(p);
                (*dfunc)(p,xi);
                dgg=gg=0.0;
                for (j=1;j<=n;j++) {
                        gg += g[j]*g[j];
                        dgg += (xi[j]+g[j])*xi[j];
                }
                if (gg == 0.0) {
                        FREEALL
                        return;
                }
                gam=dgg/gg;
                for (j=1;j<=n;j++) {
                        g[j] = -xi[j];
                        xi[j]=h[j]=g[j]+gam*h[j];
                }
        }
        nrerror("Too many iterations in frprmn");
}
#undef ITMAX
#undef EPS
#undef FREEALL
#undef NRANSI
// (C) Copr. 19 

// ## End of functions for the optimization


// ## Functions from nrutil.c

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>

#define NR_END 1
#define FREE_ARG char*

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
        fprintf(stderr,"Numerical Recipes run-time error...\n");
        fprintf(stderr,"%s\n",error_text);
        fprintf(stderr,"...now exiting to system...\n");
        exit(1);
}



float *nrvector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
        float *v;

        v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
        if (!v) nrerror("allocation failure in nrvector()");
        return v-nl+NR_END;
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
        int *v;

        v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
        if (!v) nrerror("allocation failure in ivector()");
        return v-nl+NR_END;
}

unsigned char *cvector(long nl, long nh)
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
        unsigned char *v;

        v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
        if (!v) nrerror("allocation failure in cvector()");
        return v-nl+NR_END;
}

unsigned long *lvector(long nl, long nh)
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
        unsigned long *v;

        v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
        if (!v) nrerror("allocation failure in lvector()");
        return v-nl+NR_END;
}

double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
        double *v;

        v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
        if (!v) nrerror("allocation failure in dvector()");
        return v-nl+NR_END;
}

float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
        float **m;

        /* allocate pointers to rows */
        m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m += NR_END;
        m -= nrl;

        /* allocate rows and set pointers to them */
        m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
        if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
        m[nrl] += NR_END;
        m[nrl] -= ncl;

        for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

        /* return pointer to array of pointers to rows */
        return m;
}

double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
        double **m;

        /* allocate pointers to rows */
        m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m += NR_END;
        m -= nrl;

        /* allocate rows and set pointers to them */
        m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
        if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
        m[nrl] += NR_END;
        m[nrl] -= ncl;

        for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

        /* return pointer to array of pointers to rows */
        return m;
}


cplx **cmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
        cplx **m;

        /* allocate pointers to rows */
        m=(cplx **) malloc((size_t)((nrow+NR_END)*sizeof(cplx*)));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m += NR_END;
        m -= nrl;

        /* allocate rows and set pointers to them */
        m[nrl]=(cplx *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(cplx)));
        if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
        m[nrl] += NR_END;
        m[nrl] -= ncl;

        for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

        /* return pointer to array of pointers to rows */
        return m;
}

int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
        int **m;

        /* allocate pointers to rows */
        m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m += NR_END;
        m -= nrl;


        /* allocate rows and set pointers to them */
        m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
        if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
        m[nrl] += NR_END;
        m[nrl] -= ncl;

        for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

        /* return pointer to array of pointers to rows */
        return m;
}

float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
        long newrl, long newcl)
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
        long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
        float **m;

        /* allocate array of pointers to rows */
        m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
        if (!m) nrerror("allocation failure in submatrix()");
        m += NR_END;
        m -= newrl;

        /* set pointers to rows */
        for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;

        /* return pointer to array of pointers to rows */
        return m;
}

float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
        long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
        float **m;

        /* allocate pointers to rows */
        m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
        if (!m) nrerror("allocation failure in convert_matrix()");
        m += NR_END;
        m -= nrl;

        /* set pointers to rows */
        m[nrl]=a-ncl;
        for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
        /* return pointer to array of pointers to rows */
        return m;
}

float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
        long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
        float ***t;

        /* allocate pointers to pointers to rows */
        t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
        if (!t) nrerror("allocation failure 1 in f3tensor()");
        t += NR_END;
        t -= nrl;

        /* allocate pointers to rows and set pointers to them */
        t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
        if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
        t[nrl] += NR_END;
        t[nrl] -= ncl;

        /* allocate rows and set pointers to them */
        t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
        if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
        t[nrl][ncl] += NR_END;
        t[nrl][ncl] -= ndl;

        for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
        for(i=nrl+1;i<=nrh;i++) {
                t[i]=t[i-1]+ncol;
                t[i][ncl]=t[i-1][ncl]+ncol*ndep;
                for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
        }

        /* return pointer to array of pointers to rows */
        return t;
}


void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
        free((FREE_ARG) (v+nl-NR_END));
}

void free_nrvector(float *v, long nl, long nh)
/* free a float nrvector allocated with nrvector() */
{
        free((FREE_ARG) (v+nl-NR_END));
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
        free((FREE_ARG) (v+nl-NR_END));
}

void free_cvector(unsigned char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
        free((FREE_ARG) (v+nl-NR_END));
}

void free_lvector(unsigned long *v, long nl, long nh)
/* free an unsigned long vector allocated with lvector() */
{
        free((FREE_ARG) (v+nl-NR_END));
}


void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float matrix allocated by matrix() */
{
        free((FREE_ARG) (m[nrl]+ncl-NR_END));
        free((FREE_ARG) (m+nrl-NR_END));
}

void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
        free((FREE_ARG) (m[nrl]+ncl-NR_END));
        free((FREE_ARG) (m+nrl-NR_END));
}


void free_cmatrix(cplx **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
        free((FREE_ARG) (m[nrl]+ncl-NR_END));
        free((FREE_ARG) (m+nrl-NR_END));
}

void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an int matrix allocated by imatrix() */
{
        free((FREE_ARG) (m[nrl]+ncl-NR_END));
        free((FREE_ARG) (m+nrl-NR_END));
}

void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a submatrix allocated by submatrix() */
{
        free((FREE_ARG) (b+nrl-NR_END));
}

void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a matrix allocated by convert_matrix() */
{
        free((FREE_ARG) (b+nrl-NR_END));
}

void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
        long ndl, long ndh)
/* free a float f3tensor allocated by f3tensor() */
{
        free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
        free((FREE_ARG) (t[nrl]+ncl-NR_END));
        free((FREE_ARG) (t+nrl-NR_END));
}






#include <math.h>
#define NRANSI
#include "nrutil.h"
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

// ## Changed by Ivan in 01/2001 from float to double

void gaussj(double **a, int n, double **b, int m) //
{
        int *indxc,*indxr,*ipiv;
        int i,icol,irow,j,k,l,ll;
        double big,dum,pivinv,temp;  //

        indxc=ivector(1,n);
        indxr=ivector(1,n);
        ipiv=ivector(1,n);
        for (j=1;j<=n;j++) ipiv[j]=0;
        for (i=1;i<=n;i++) {
                big=0.0;
                for (j=1;j<=n;j++)
                        if (ipiv[j] != 1)
                                for (k=1;k<=n;k++) {
                                        if (ipiv[k] == 0) {
                                                if (abs(a[j][k]) >= big) {
                                                        big=abs(a[j][k]);
                                                        irow=j;
                                                        icol=k;
                                                }
                                        } else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
                                }
                ++(ipiv[icol]);
                if (irow != icol) {
                        for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
                        for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
                }
                indxr[i]=irow;
                indxc[i]=icol;
                if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
                pivinv=1.0/a[icol][icol];
                a[icol][icol]=1.0;
                for (l=1;l<=n;l++) a[icol][l] *= pivinv;
                for (l=1;l<=m;l++) b[icol][l] *= pivinv;
                for (ll=1;ll<=n;ll++)
                        if (ll != icol) {
                                dum=a[ll][icol];
                                a[ll][icol]=0.0;
                                for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
                                for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
                        }
        }
        for (l=n;l>=1;l--) {
                if (indxr[l] != indxc[l])
                        for (k=1;k<=n;k++)
                                SWAP(a[k][indxr[l]],a[k][indxc[l]]);
        }
        free_ivector(ipiv,1,n);
        free_ivector(indxr,1,n);
        free_ivector(indxc,1,n);
}
#undef SWAP
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). */






#include <math.h>
#define NRANSI
#include "nrutil.h"
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

// ## Changed by Ivan in 04/2001 from double to cplx
void gaussj_complex(cplx **a, int n, cplx **b, int m) //
{
        int *indxc,*indxr,*ipiv;
        int i,icol,irow,j,k,l,ll;
        cplx dum,pivinv,temp;  //
        double big;
        indxc=ivector(1,n);
        indxr=ivector(1,n);
        ipiv=ivector(1,n);
        for (j=1;j<=n;j++) ipiv[j]=0;
        for (i=1;i<=n;i++) {
                big=0.0;
                for (j=1;j<=n;j++)
                        if (ipiv[j] != 1)
                                for (k=1;k<=n;k++) {
                                        if (ipiv[k] == 0) {
                                                if (abs(a[j][k]) >= big) {
                                                        big=abs(a[j][k]);
                                                        irow=j;
                                                        icol=k;
                                                }
                                        } else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
                                }
                ++(ipiv[icol]);
                if (irow != icol) {
                        for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
                        for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
                }
                indxr[i]=irow;
                indxc[i]=icol;
                if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
                pivinv=1.0/a[icol][icol];
                a[icol][icol]=1.0;
                for (l=1;l<=n;l++) a[icol][l] *= pivinv;
                for (l=1;l<=m;l++) b[icol][l] *= pivinv;
                for (ll=1;ll<=n;ll++)
                        if (ll != icol) {
                                dum=a[ll][icol];
                                a[ll][icol]=0.0;
                                for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
                                for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
                        }
        }
        for (l=n;l>=1;l--) {
                if (indxr[l] != indxc[l])
                        for (k=1;k<=n;k++)
                                SWAP(a[k][indxr[l]],a[k][indxc[l]]);
        }
        free_ivector(ipiv,1,n);
        free_ivector(indxr,1,n);
        free_ivector(indxc,1,n);
}
#undef SWAP
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). */


float erffc(float x)
{
        float gammp(float a, float x);
        float gammq(float a, float x);

        return x < 0.0 ? 1.0+gammp(0.5,x*x) : gammq(0.5,x*x);
}
/* (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). */

float gammp(float a, float x)
{
        void gcf(float *gammcf, float a, float x, float *gln);
        void gser(float *gamser, float a, float x, float *gln);
        void nrerror(char error_text[]);
        float gamser,gammcf,gln;

        if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp");
        if (x < (a+1.0)) {
                gser(&gamser,a,x,&gln);
                return gamser;
        } else {
                gcf(&gammcf,a,x,&gln);
                return 1.0-gammcf;
        }
}
/* (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). */

#include <math.h>
#define ITMAX 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30

void gcf(float *gammcf, float a, float x, float *gln)
{
        float gammln(float xx);
        void nrerror(char error_text[]);
        int i;
        float an,b,c,d,del,h;

        *gln=gammln(a);
        b=x+1.0-a;
        c=1.0/FPMIN;
        d=1.0/b;
        h=d;
        for (i=1;i<=ITMAX;i++) {
                an = -i*(i-a);
                b += 2.0;
                d=an*d+b;
                if (fabs(d) < FPMIN) d=FPMIN;
                c=b+an/c;
                if (fabs(c) < FPMIN) c=FPMIN;
                d=1.0/d;
                del=d*c;
                h *= del;
                if (fabs(del-1.0) < EPS) break;
        }
        if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");
        *gammcf=exp(-x+a*log(x)-(*gln))*h;
}
#undef ITMAX
#undef EPS
#undef FPMIN
/* (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). */

#include <math.h>
#define ITMAX 100
#define EPS 3.0e-7

void gser(float *gamser, float a, float x, float *gln)
{
        float gammln(float xx);
        void nrerror(char error_text[]);
        int n;
        float sum,del,ap;

        *gln=gammln(a);
        if (x <= 0.0) {
                if (x < 0.0) nrerror("x less than 0 in routine gser");
                *gamser=0.0;
                return;
        } else {
                ap=a;
                del=sum=1.0/a;
                for (n=1;n<=ITMAX;n++) {
                        ++ap;
                        del *= x/ap;
                        sum += del;
                        if (fabs(del) < fabs(sum)*EPS) {
                                *gamser=sum*exp(-x+a*log(x)-(*gln));
                                return;
                        }
                }
                nrerror("a too large, ITMAX too small in routine gser");
                return;
        }
}
#undef ITMAX
#undef EPS
/* (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). */


#include <math.h>

float gammln(float xx)
{
        double x,y,tmp,ser;
        static double cof[6]={76.18009172947146,-86.50532032941677,
                24.01409824083091,-1.231739572450155,
                0.1208650973866179e-2,-0.5395239384953e-5};
        int j;

        y=x=xx;
        tmp=x+5.5;
        tmp -= (x+0.5)*log(tmp);
        ser=1.000000000190015;
        for (j=0;j<=5;j++) ser += cof[j]/++y;
        return -tmp+log(2.5066282746310005*ser/x);
}
/* (C) Copr. 1986-92 Numeric*/


float gammq(float a, float x)
{
        void gcf(float *gammcf, float a, float x, float *gln);
        void gser(float *gamser, float a, float x, float *gln);
        void nrerror(char error_text[]);
        float gamser,gammcf,gln;

        if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq");
        if (x < (a+1.0)) {
                gser(&gamser,a,x,&gln);
                return 1.0-gamser;
        } else {
                gcf(&gammcf,a,x,&gln);
                return gammcf;
        }
}
/* (C) Copr. 1986-92 Numerical Recipes Software *1%:?!$:V3+). */



// ##### End of functions obtained from Numerical Recipes

---