laser_field.cc

#include <mpi.h>
#include <cmath>
#include "laser_field.h"
#include "g_index.h"
#include "g_gas.h"
#include "conf.h"
#include "Paralleler/paral_parent.h"
#include "g_mesh.h"
#include "g_flow.h"
#include "g_unsteady.h"
#include "g_eint.h"
#include "flowf.h"

extern cgUnsteady * Unsteady;
extern cgIndex  * Index;    // Indexation
extern cgGases  *Gases;   // Gases list !new
extern cConf * Conf; // Calculation configuration
extern ParalParent * paralleler; // Paralleler father
extern cgMesh * Mesh;   // Calculation mesh
extern cgEInt *EIntMod; // model of internal energy
extern cgFlowField *Field;	// FlowField

Laser * LaserL = NULL;		//Left Laser
Laser * LaserR = NULL;		//Right Laser
Laser_Field * LField = NULL;	//Electrical laser field
MacroParameters * MParams = NULL;	//Macroparameters Class object

const R8 CSpeed = 3.0e8;

template <typename T>
    std::string to_string(T value)
    {
      //create an output string stream
      std::ostringstream os ;

      //throw the value into the string stream
      os << value ;

      //convert the string stream into a string and return
      return os.str() ;
    }

Laser::Laser()
{
	LaserEnergy = 0.0;			read.Add("Energy", &LaserEnergy);
	LaserWL = 0.0;				read.Add("WL", &LaserWL);
	LaserDiam = 0.0;			read.Add("Diam", &LaserDiam);
	LaserPulseWidth = 0.0;		read.Add("PulseWidth", &LaserPulseWidth);
	LaserPulseCenter = 0.0;		read.Add("PulseCenter", &LaserPulseCenter);
	LaserYPosition = 0.0;		read.Add("YPosition", &LaserYPosition);
	LaserZPosition = 0.0;		read.Add("ZPosition", &LaserZPosition);
}

void Laser::Make()
{
	//call recomputation or additional info if necessary
	Intensity = 0.83*LaserEnergy/(LaserDiam*LaserDiam)/LaserPulseWidth;		//doesn't include time and radial exp() dependencies 
	LaserK = 2.0*3.1415926/LaserWL;
	LaserOmega = LaserK*CSpeed;
	if (paralleler->GetIProc() == 0) printf("Laser wl is: %le\n", LaserWL);
	if (paralleler->GetIProc() == 0) printf("Laser Omega is: %le\n", LaserOmega);
	if (paralleler->GetIProc() == 0) printf("Laser Intensity is: %le\n", Intensity);

	return;
}

void Laser_Field::Make()
{
	BulkLength = Mesh->dcell * sqrt(Mesh->Vol(0)); //length if cells are the same
	if (paralleler->GetIProc() == 0) printf("Bulk length is: %le [m]\n", BulkLength);

	ForceX = new R8[Mesh->dcell];
	for(int i = 0; i < Mesh->dcell; i++) ForceX[i] = 0.0;

	FullEField = new R8[Mesh->dcell];
	for(int i = 0; i < Mesh->dcell; i++) FullEField[i] = 0.0;

	EL_plus = new complex<double>[Mesh->dcell];
	EL_minus = new complex<double>[Mesh->dcell]; 
	ER_plus = new complex<double>[Mesh->dcell];
  	ER_minus = new complex<double>[Mesh->dcell];
  	E_Interference = new double[Mesh->dcell];

	if (paralleler->GetIProc() == 1) file = fopen("file_laser.txt","w");

	if (paralleler->GetIProc() == 0) cout << "Lattice Speed is " << (LaserL->LaserOmega - LaserR->LaserOmega)/(LaserL->LaserK + LaserR->LaserK) << endl;

	count_field_write = 0;
	return;
}

R8 Laser_Field::ExternalForceX(R8 pos) //goes to VelocityChange
{	
	int NumForceXCell = int(pos/BulkLength * Mesh->dcell);
	
	if (NumForceXCell > int(Mesh->dcell)) NumForceXCell = Mesh->dcell;
	if (NumForceXCell < 0) NumForceXCell = 0;

	R8 local_ForceX = ForceX[NumForceXCell];

	//if Intensity is time dependent
	//R8 curTime = Conf->time - Conf->Step;
	//local_ForceX = local_ForceX*exp(-2.77*pow((curTime-LaserL->LaserPulseCenter)/LaserL->LaserPulseWidth,2));
	
	return local_ForceX;
}

void Laser_Field::VelocityChange(V8 * vel, V8 * pos, R8 Step, int gas, int cell) //works from move2d_rect_pb.cc
{
	R8 Gmass = Gases->Get(gas)->mass;	//get gas particle mass
	
	//Euler method (old) v + F/m * dt
	//vel->x = vel->x + ExternalForceX(pos->x)/Gmass*Step;

	//RK-4 method (new) 
	R8 k0 = Step*vel->x;
	R8 q0 = Step*ExternalForceX(pos->x)/Gmass;

	R8 k1 = Step*(vel->x + q0/2.0);
	R8 q1 = Step*ExternalForceX(pos->x+k0/2.0)/Gmass;

	R8 k2 = Step*(vel->x + q1/2.0);
	R8 q2 = Step*ExternalForceX(pos->x+k1/2.0)/Gmass;	
	
	R8 k3 = Step*(vel->x + q2);
	R8 q3 = Step*ExternalForceX(pos->x+k2)/Gmass;

	vel->x = vel->x + (q0 + 2.0*q1 + 2.0*q2 + q3)/6.0;

	return;
}

void Laser_Field::ParallelCollect() 	//called from dsmc.cc
{

	MParams->DenstGather();		//Gathering density on all procs
	MParams->refractIndexGather(); //calc refr index on all procs
	return;
}

void Laser_Field::CalcForce_Analytical()
{
	MParams->GetDenstForm(); // !moved from DenstGather

	Sparam = abs(MParams->DenAmplitude)*molAlpha/(2.0*epsilon0)*sqrt(LaserR->LaserOmega*LaserL->LaserOmega)/(2*CSpeed); //see Shneider OPTICS 17a, 17b
	//if (paralleler->GetIProc() == 1) printf("Sparam: %le\n", Sparam);

	R8 curTime = Conf->time - Conf->Step;					//current time on the beggining of new iteration	
	R8 q = (LaserR->LaserK + LaserL->LaserK);				//k1 + k2
	R8 omega = LaserR->LaserOmega - LaserL->LaserOmega;		//w1 - w2
	R8 pos;			// pos equals average x position in current cell

	R8 Int1 = LaserL->Intensity;	//laserR intensity
	R8 Int2 = LaserR->Intensity; 	//laserL intensity
	R8 local_ForceX;

	for(int i = 0; i < Mesh->dcell; i++) ForceX[i] = 0;	

	for(int i = 0; i < paralleler->cells_number_on_processor(); i++)
	{
		pos = BulkLength / Mesh->dcell * (paralleler->global_number_from_local(i,paralleler->GetIProc()) + 1/2);
				
		Int1 = 1.0/cosh(Sparam*BulkLength) * sqrt((pow(LaserL->Intensity*cosh(Sparam*(pos - BulkLength)),2.0) + pow(LaserR->Intensity*sinh(Sparam*pos),2.0)*pow(LaserR->LaserOmega/LaserL->LaserOmega,3.0))); //plus analytical nonlinear
		Int2 = 1.0/cosh(Sparam*BulkLength) * sqrt((pow(LaserR->Intensity*cosh(Sparam*pos),2.0) + pow(LaserL->Intensity*sinh(Sparam*(pos - BulkLength)),2.0)*pow(LaserL->LaserOmega/LaserR->LaserOmega,3.0))); //plus analytical nonlinear

		local_ForceX = (-molAlpha/CSpeed/epsilon0*q) * sqrt(Int1 * Int2) * sin(q*pos-omega*curTime);
		//if Intensity is time dependent
		//local_ForceX = local_ForceX*exp(-2.77*pow((curTime-LaserL->LaserPulseCenter)/LaserL->LaserPulseWidth,2));

		ForceX[paralleler->global_number_from_local(i,paralleler->GetIProc())] = local_ForceX;
	}	

	paralleler->sum_by_all_processors(ForceX, Mesh->dcell);	//sums force array over the all procs
}
//time averaged interference of two counterpropogating fields
void Laser_Field::CalcInterference() //called from CalcField_LayerModel
{
	R8 curTime = Conf->time - Conf->Step;	//current time on the beggining of new iteration	
	complex<double> timePhase0 (0.0, LaserL->LaserOmega * curTime);
	complex<double> timePhase1 (0.0, LaserR->LaserOmega * curTime);
	int size = Mesh->dcell;
	double dx = BulkLength / Mesh->dcell;
	double k0 = LaserL->LaserK;
	double k1 = LaserR->LaserK;


	for (int i = 0; i < Mesh->dcell; i++){
		complex<double> spatPhase0 (0.0, *(MParams->refractIndex+i)*k0*(double(i)-double(size-1)/2.0)*dx);
		complex<double> spatPhase1 (0.0, *(MParams->refractIndex+i)*k1*(double(i)-double(size-1)/2.0)*dx);

		*(E_Interference+i) = 0.5 * pow(abs( *(EL_plus+i)*exp(spatPhase0-timePhase0) + *(ER_plus+i) * exp(spatPhase0-timePhase0) + *(EL_minus+i)*exp(-spatPhase1-timePhase1)  + *(ER_minus+i)*exp(-spatPhase1-timePhase1) ) , 2.0);
	}//!

	
	return;
}

void Laser_Field::CalcField_LayerModel() //called from CalcForce_LayerModel
{
	int size = Mesh->dcell;

	vector<cgUnsteadyRegion *>::iterator it;
 	it =Unsteady->regions.begin();
	double DensInitial = (*it)->flows->ndens;
	
	//1. Calc right, left fields and normalyze them (based on initial intensities)
	complex<double> Eleft_plus;
	complex<double> Eleft_minus;
	complex<double> Eright_plus;
	complex<double> Eright_minus;
	
	double nr, nl;
	double nr_externBoundary = sqrt(1.0+DensInitial*molAlpha/epsilon0);
	double nl_externBoundary = sqrt(1.0+DensInitial*molAlpha/epsilon0);
	double k0 = LaserL->LaserK;
	double k1 = LaserR->LaserK;
	double dx = BulkLength / Mesh->dcell;

	
	//==============================
	//Starting left wave calculation
	Eright_plus = complex<double>(1.0,0.0);		//Волна, которая вылетает справа из среды
	Eright_minus = complex<double>(0.0,0.0);	//Волна, которая влетает в среду справа (нулевая)

	for (int i = size; i >= 0; i--)
	{
		if (i == size)
			nr = nr_externBoundary;
		else
			nr = *(MParams->refractIndex + i);

		if (i == 0)
			nl = nl_externBoundary;
		else
			nl = *(MParams->refractIndex + i - 1);

		complex<double> phaseR (0.0, k0*nr*(double(i)-double(size)/2.0)*dx); // i * kR * z
		complex<double> phaseL (0.0, k0*nl*(double(i)-double(size)/2.0)*dx); // i * kL * z

		Eleft_plus = (Eright_plus*exp(phaseR)*(nr+nl)/(2.0*nl) + Eright_minus*exp(-phaseR)*(nl-nr)/(2.0*nl)) * exp(-phaseL); 
		Eleft_minus = (Eright_plus*exp(phaseR)*(nl-nr)/(2.0*nl) + Eright_minus*exp(-phaseR)*(nr+nl)/(2.0*nl)) * exp(phaseL);

		
		//заполнение и сохранение массивов
		if (i > 0){
			*(EL_plus+i-1) = Eleft_plus;
			*(EL_minus+i-1) = Eleft_minus;
			}
		

		//переход на след границу
		Eright_plus = Eleft_plus;
		Eright_minus = Eleft_minus;
	}
	//Normalization of the left wave (so E+[-inf] = (E0,0) intensity based amplitude)
	for (int i = size-1; i >= 0; i--)
	{	
		*(EL_minus+i) = *(EL_minus+i) / Eleft_plus * sqrt(2.0*LaserL->Intensity / sqrt(epsilon0/mu0 * (1.0 + molAlpha*DensInitial/epsilon0)));
		*(EL_plus+i) = *(EL_plus+i) / Eleft_plus  * sqrt(2.0*LaserL->Intensity / sqrt(epsilon0/mu0 * (1.0 + molAlpha*DensInitial/epsilon0)));
	}
	//Left wave calculation completed!

	//==============================
	//Starting right wave calculation
	Eleft_plus = complex<double>(0.0,0.0);	//Волна, которая вылетает справа из среды 
	Eleft_minus = complex<double>(1.0,0.0);	//Волна, которая влетает в среду справа (нулевая)

	for (int i = 0; i <= size; i++)
	{
		if (i == 0)
			nl = nl_externBoundary;
		else
			nl = *(MParams->refractIndex + i - 1);

		if (i == size)
			nr = nr_externBoundary;
		else
			nr = *(MParams->refractIndex + i);
		
		complex<double> phaseR (0.0, k0*nr*(double(i)-double(size)/2.0)*dx); // i * kR * z
		complex<double> phaseL (0.0, k0*nl*(double(i)-double(size)/2.0)*dx); // i * kL * z

		Eright_plus = (Eleft_plus*exp(phaseL)*(nr+nl)/(2.0*nr) + Eleft_minus*exp(-phaseL)*(nr-nl)/(2.0*nr)) * exp(-phaseR);
		Eright_minus = (Eleft_plus*exp(phaseL)*(nr-nl)/(2.0*nr) + Eleft_minus*exp(-phaseL)*(nr+nl)/(2.0*nr)) * exp(phaseR);
	
		if (i < size){
			*(ER_plus+i) = Eright_plus;
			*(ER_minus+i) = Eright_minus;
		}
			
		Eleft_plus = Eright_plus;
		Eleft_minus = Eright_minus;
	}

	//Normalization of the right wave (so E-[inf] = (E0,0) intensity based amplitude)
	for (int i = 0; i < size; i++)
	{
		*(ER_plus+i) = *(ER_plus+i) / Eright_minus  * sqrt(2.0*LaserR->Intensity / sqrt(epsilon0/mu0 * (1.0 + molAlpha*DensInitial/epsilon0)));
		*(ER_minus+i) = *(ER_minus+i) / Eright_minus  * sqrt(2.0*LaserR->Intensity / sqrt(epsilon0/mu0 * (1.0 + molAlpha*DensInitial/epsilon0)));
	}
	//Right wave calculation completed!
	//==============================
	
	// calc <E^2> = 1/2 |E+ +E-|^2
	CalcInterference();

	//2. Calc full field U = -1/2 alpha <E^2>
	for (int i = 0; i < size; i++) *(FullEField+i) = *(E_Interference+i) * (-0.5) * molAlpha; 
	
	return;
}

void Laser_Field::CalcForce_LayerModel()//called from CalcForce
{
	//First we need to calculate full laser field based on a layer model
	for(int i = 0; i < Mesh->dcell; i++) ForceX[i] = 0;

	CalcField_LayerModel();

	for(int i = 0; i < paralleler->cells_number_on_processor(); i++)
	{
		int i_global_ = paralleler->global_number_from_local(i,paralleler->GetIProc());
		
		if (i_global_ == 0){ 
			double x1 = double(i_global_) * BulkLength / Mesh->dcell;
			double x2 = double(i_global_ + 1) * BulkLength / Mesh->dcell;
			double x3 = double(i_global_ + 2) * BulkLength / Mesh->dcell;
			double y1 = FullEField[i_global_];
			double y2 = FullEField[i_global_ + 1];
			double y3 = FullEField[i_global_ + 2];
			double a = (y3 - (x3*(y2-y1)+x2*y1-x1*y2)/(x2-x1) )/(x3*(x3 - x2 - x1) + x1*x2);
			double b = (y2-y1)/(x2-x1) - (x1+x2)*a;
			ForceX[i_global_] = -(2*a*x1 + b);
			}
		else if (i_global_ == Mesh->dcell - 1){
			double x1 = double(i_global_-2) * BulkLength / Mesh->dcell;
			double x2 = double(i_global_ - 1) * BulkLength / Mesh->dcell;
			double x3 = double(i_global_) * BulkLength / Mesh->dcell;
			double y1 = FullEField[i_global_-2];
			double y2 = FullEField[i_global_ - 1];
			double y3 = FullEField[i_global_];
			double a = (y3 - (x3*(y2-y1)+x2*y1-x1*y2)/(x2-x1) )/(x3*(x3 - x2 - x1) + x1*x2);
			double b = (y2-y1)/(x2-x1) - (x1+x2)*a;
			ForceX[i_global_] = -(2*a*x3 + b);
			}  
		else  //central approxim F'(x)
			ForceX[i_global_] = -(FullEField[i_global_+1] - FullEField[i_global_-1]) / (2.0 * BulkLength / Mesh->dcell);
	}

	paralleler->sum_by_all_processors(ForceX, Mesh->dcell);	//sums force array over the all procs

	return;
}
// after paralleler calculate force gradients
void Laser_Field::CalcForce() // called from dsmc.cc
{
	AnalytModel = false;
	LayerModel = true;
	
	if (AnalytModel == true)
		CalcForce_Analytical(); 
	else if (LayerModel == true)
		CalcForce_LayerModel();

	if (paralleler->GetIProc() == 0){
		count_field_write++;
		if (count_field_write % 10 == 0){

		std::string fileName;
		fileName = "Field/field_" + to_string( Conf->time / Conf->Step + 0.0 ) + ".txt";
		const char *cstr = fileName.c_str();
		FILE * fileE = fopen(cstr, "w");
		//UNCOMMENT FOR FIELD
		//for (int i = 0; i < Mesh->dcell; i++) fprintf(fileE, "%e %e %e %e %e\n", *(MParams->denst+i), *(ForceX+i), *(E_Interference+i), abs (*(EL_plus+i) + (*(ER_plus+i))), abs( *(EL_minus+i) + (*(ER_minus+i))) );
		for (int i = 0; i < Mesh->dcell; i++) fprintf(fileE, "%e %e %e %e %e\n", *(MParams->denst+i), abs (*(EL_minus+i)), abs (*(EL_plus+i)), abs (*(EL_plus+i) + (*(ER_plus+i))),  abs( *(EL_minus+i) + (*(ER_minus+i))) );
		fclose(fileE);
		}
		}
	
	if (paralleler->GetIProc() == 2) { MParams->GetDenstForm( int(Mesh->dcell * 3 / 5 ) );}
	if (paralleler->GetIProc() == 1) { MParams->GetDenstForm( int(Mesh->dcell * 2 / 5 ) );}
	if (paralleler->GetIProc() == 0) {MParams->GetDenstForm( int(Mesh->dcell / 2 ) ); }
	

	double MeanTempX = 0.0;
	double MeanTempY = 0.0;
	double MeanTempZ = 0.0;
	double MeanDensity = 0.0;
	double MeanUx = 0.0;
	double MeanUy = 0.0;
	double MeanUz = 0.0;
	double MeanT_Temp = 0.0;
	//for(int i = 0; i < paralleler->cells_number_on_processor(); i++) MeanTemp = MeanTemp + EIntMod->GetMeanTTrn( i );
	//paralleler->sum_by_all_processors(&MeanTemp, sizeof(double));
	//MeanTemp = MeanTemp / Mesh->dcell;
	for (int i = int(Mesh->dcell*0/20); i < int(Mesh->dcell*20/20); i++)
	{
		if (paralleler->whose_is_cell(i) == paralleler->GetIProc()){
			double Ux = Field->GetValue("Ux", paralleler->local_number_from_global(i), -1, 0);
			double Uy = Field->GetValue("Uy", paralleler->local_number_from_global(i), -1, 0);
			double Uz = Field->GetValue("Uz", paralleler->local_number_from_global(i), -1, 0);

			double Tx = Field->GetValue("Tx", paralleler->local_number_from_global(i), -1, 0);
			double Ty = Field->GetValue("Ty",  paralleler->local_number_from_global(i), -1, 0);
			double Tz = Field->GetValue("Tz",  paralleler->local_number_from_global(i), -1, 0);			

			MeanTempX = MeanTempX + (Tx  + 6.63e-26 * Ux*Ux/3.0/1.38e-23       )*  MParams->denst[i];
			MeanTempY = MeanTempY + (Ty  + 6.63e-26 * Uy*Uy/3.0/1.38e-23       )*  MParams->denst[i];
			MeanTempZ = MeanTempZ + (Tz  + 6.63e-26 * Uz*Uz/3.0/1.38e-23       )*  MParams->denst[i];
			MeanDensity = MeanDensity + MParams->denst[i];
			MeanUx = MeanUx + Ux *  MParams->denst[i]; 
			MeanUy = MeanUy + Uy *  MParams->denst[i]; 
			MeanUz = MeanUz + Uz *  MParams->denst[i]; 

			MeanT_Temp = MeanT_Temp + (Tx + Ty + Tz)/3.0 * MParams->denst[i];
			}
	}
	paralleler->sum(MeanTempX);
	paralleler->sum(MeanTempY);
	paralleler->sum(MeanTempZ);	
	paralleler->sum(MeanDensity);
	paralleler->sum(MeanUx);
	paralleler->sum(MeanUy);
	paralleler->sum(MeanUz);

	paralleler->sum(MeanT_Temp);

	MeanUx = MeanUx / MeanDensity;
	MeanUy = MeanUy / MeanDensity;
	MeanUz = MeanUz / MeanDensity;
	MeanTempX = MeanTempX/ MeanDensity - 6.63e-26*MeanUx*MeanUx/3.0/1.38e-23;
	MeanTempY = MeanTempY/ MeanDensity - 6.63e-26*MeanUy*MeanUy/3.0/1.38e-23;
	MeanTempZ = MeanTempZ/ MeanDensity - 6.63e-26*MeanUz*MeanUz/3.0/1.38e-23;
	
	MeanT_Temp = MeanT_Temp/ MeanDensity;

	if (paralleler->GetIProc() == 0) {fprintf(MParams->fileT, "%le %le %le\n", MeanT_Temp, MeanTempX,  MeanUx);}
	
	return;
}

void MacroParameters::Make()
{
	denst = new R8[Mesh->dcell];	//creates new array for denst sampling over the whole area
	refractIndex = new R8[Mesh->dcell]; //creates new array for refractive index values

	EffWL = LaserL->LaserWL*LaserR->LaserWL/(LaserL->LaserWL+LaserR->LaserWL);	//effective wavelength wl = wl1*wl2 / (wl1+wl2)
	EffK = LaserL->LaserK + LaserR->LaserK;
	EffOmega = LaserL->LaserOmega - LaserR->LaserOmega; 
	//printf("%le\n",EffOmega/EffK);

	ArrSize = 10*int(EffWL/sqrt(Mesh->Vol(0)));	//number of cells for sin-like density calculation
	//printf("%d\n", ArrSize);

	DenPhase = 3.1415926/2.0;
	if (paralleler->GetIProc() == 0) {
		std::string fileName;
		fileName = "ampl_phase_" + to_string(paralleler->GetIProc() ) + ".txt";
		const char *cstr = fileName.c_str();
		file = fopen(cstr,"w");

		fileT = fopen("temperature.txt", "w");
	}
	if (paralleler->GetIProc() == 1) {
		std::string fileName;
		fileName = "ampl_phase_" + to_string(paralleler->GetIProc() ) + ".txt";
		const char *cstr = fileName.c_str();
		file = fopen(cstr,"w");
	}
	if (paralleler->GetIProc() == 2) {
		std::string fileName;
		fileName = "ampl_phase_" + to_string(paralleler->GetIProc() ) + ".txt";
		const char *cstr = fileName.c_str();
		file = fopen(cstr,"w");
	}


	return;
}

void MacroParameters::refractIndexGather() //called from Laser_Field::ParallelCollect
{
	for (int i = 0; i < Mesh->dcell; i++) refractIndex[i] = 1.0;

	//n(x) = sqrt(1 + N(x)*alpha/eps0)
	for (int i = 0; i < Mesh->dcell; i++) refractIndex[i] = sqrt(1.0 + denst[i] * LField->molAlpha / LField->epsilon0);	
	
	return;
}

void MacroParameters::DenstGather()
{
	int cells_number = paralleler->cells_number_on_processor();	//number of cells on current proc
	R8 NumPtc = 0;

	int curProc = paralleler->GetIProc();
	int global_i;

	for(int i = 0; i < Mesh->dcell; i++) denst[i] = 0;	

	for(int i = 0; i < cells_number; i++)
	{
		NumPtc = Index->NumPtc(i);	//i - local number, returns number of particles
		global_i = paralleler->global_number_from_local(i, curProc);
		denst[global_i] = NumPtc*Conf->FNum/Mesh->Vol(i);	//number of particles * FNum / CellVolume
	}
	paralleler->sum_by_all_processors(denst, Mesh->dcell);	//sums denst array over the all procs
	return;
}

void MacroParameters::GetDenstForm() //to optimize
{
	GetDenstForm(0);
	return;
}

void MacroParameters::GetDenstForm(int CellNumber) //to optimize
{

	StartCell = CellNumber;

	R8 phase;	//for phase calculation
	R8 phase2 = 0.0;
	R8 curPhase = 0.0;
	R8 phaseStep = 3.1415926*sqrt(2.0)/100.0;
	phase = DenPhase - sqrt(2.0)*phaseStep;

	R8 DenstAverage = 0.0;

	//StartCell = int(Mesh->dcell / 2 );//0
	//if (paralleler->GetIProc() == 1) printf("%d\n",StartCell);

	/*Calculating average number density */
	for (int i = StartCell; i < StartCell+ArrSize; i++) DenstAverage += denst[i];
	DenAverage = DenstAverage/ArrSize;

	R8 error1 = 1.0;
	R8 error2 = 1.0;

	while (error1*error2 >= 0.0)
	{
		error1 = DenstDiffToMinimize(phase);
		error2 = DenstDiffToMinimize(phase+=phaseStep);
	}

	phase2 = phase;
	phase -= phaseStep;

	R8 error = DenstFromSinDeviation(phase);
	R8 errorPrev;
	int j = 0;

	do
	{
		errorPrev = error;
		curPhase = (phase2+phase)/2.0;
		phaseStep = phaseStep/2.0;

		if (DenstDiffToMinimize(phase)*DenstDiffToMinimize(curPhase) <= 0)
			{phase2 = phase;
			phase = curPhase;}
		else
			phase = curPhase;
		error = DenstFromSinDeviation(phase);
		j++;
	}
	while (abs((error - errorPrev)/error) > 1.0e-5);

	while (phase > 3.1415926) phase -= 3.1415926;

	DenPhase = phase;
	DenAmplitude = GetDenstAmplitude(phase);

	//if (paralleler->GetIProc() == 1) printf("%le\n",DenPhase);
	//if (paralleler->GetIProc() == 1) printf("%d\n",j);
	
	//if (paralleler->GetIProc() == 1)

	
	fprintf(file, "%le %le\n", DenPhase,DenAmplitude); //phaseStep

	//if (paralleler->GetIProc() == 1) printf("Sparam: %le\n", LField->Sparam);
	
	//R8 TestInt1 = 1.0/cosh(LField->Sparam*LField->BulkLength) * sqrt(pow(LaserL->Intensity*cosh(LField->Sparam*( - LField->BulkLength/2.0)),2.0) + pow(LaserR->Intensity*sinh(LField->Sparam*LField->BulkLength/2.0),2.0)*pow(LaserR->LaserOmega/LaserL->LaserOmega,3.0)); //plus analytical nonlinear
	//if (paralleler->GetIProc() == 1) fprintf(file, "%le\n", TestInt1);
	return;
}

R8 MacroParameters::DenstFromSinDeviation(R8 phase)
{
	R8 curTime = Conf->time - Conf->Step;
	R8 deviation = 0.0;
	R8 ampl = GetDenstAmplitude(phase);
	for (int i = StartCell; i < StartCell+ArrSize; i++)
		deviation += pow((ampl*sin(EffK*i*sqrt(Mesh->Vol(0)) - EffOmega*curTime + phase) - (denst[i]-DenAverage)),2.0);

	return deviation;
}

R8 MacroParameters::DenstDiffToMinimize(R8 phase)
{
	R8 curTime = Conf->time - Conf->Step;
	R8 error = 0.0;
	for (int i = StartCell; i < StartCell+ArrSize; i++)
		error += (GetDenstAmplitude(phase)*sin(EffK*i*sqrt(Mesh->Vol(0)) - EffOmega*curTime + phase) - (denst[i]-DenAverage)) * cos(EffK*i*sqrt(Mesh->Vol(0)) - EffOmega*curTime + phase);
	
	return error;
}

R8 MacroParameters::GetDenstAmplitude(R8 phase)	//calculates assumed sin-like density amplitude by the given phase
{	
	R8 curTime = Conf->time - Conf->Step;
	R8 fsin = 0.0;
	R8 sinsin = 0.0;

	for (int i = StartCell; i < StartCell+ArrSize; i++)
	{
		fsin += (denst[i]-DenAverage)*sin(EffK*i*sqrt(Mesh->Vol(0)) - EffOmega*curTime + phase);
		sinsin += pow(sin(EffK*i*sqrt(Mesh->Vol(0)) - EffOmega*curTime + phase),2.0);
	}
	return fsin/sinsin;
}

void Laser::Print(FILE * fp)
{	
	return;
}