refract_3.0.cc

#define _USE_MATH_DEFINES
#include <iostream>
#include <iomanip>
#include <cmath>
#include <stdio.h>
#include <complex>      // std::complex, std::abs
#include <cstdlib>
#include <ctime>
using namespace std;

double ampl = 3.7e-5;
double n0 = 1.0+ampl*10.0;

double PI = 3.14159265359;
double bulkLength = 4.0e2;
int discretSize = 100;
int size = int(bulkLength)*discretSize;
double latticeSpeed = 0.0;
double cspeed = 3.0e8;
double lambda0 = 1.0*double(discretSize);
double k0 = 2.0*PI/lambda0;
double lambda1 = lambda0*(1 + 2.0*latticeSpeed/cspeed*n0);
double k1 = 2.0*PI/lambda1;



void set_n(double * n, FILE * file_n)
{
	std::srand(std::time(0)); // use current time as seed for random generator
    
    cout << "lattice speed: " << cspeed*(lambda1 - lambda0) / (lambda1 + lambda0) / n0 << " m/s" << endl; 
    cout << "bulk Length is: " << bulkLength*532.0e-7 << " cm\n" << endl;
	
	/*while ( abs(cos((k0+k1)*n0*double(size)) - 1.0) > 1.0e-2)
		size++;
	if ( abs(cos((k0+k1)*n0*double(size)) - cos((k0+k1)*n0*double(size+1))) < 1.0e-2)
		size++;
	cout << "added " << size - bulkLength*discretSize << " points." << endl;*/

	int BulkCenter = size / 2;
	
	for (int i = 0; i <= BulkCenter; i++)
	{
		*(n+BulkCenter+i) = n0 + ampl*cos((k0+k1)*n0*double(i));
		*(n+BulkCenter-i) = n0 + ampl*cos((k0+k1)*n0*double(i));
	//	cout << *(n+BulkCenter+i) - *(n+BulkCenter-i) << endl;
	//	cin.get();
	}

	/*for(int i = 0; i <= size; i++){

		double random_variable = double(std::rand()) / double(RAND_MAX);

		*(n+i) = n0 + ampl*(cos((k0+k1)*n0*double(i)) + 0.0*(-2.0*random_variable+1.0));
		//*(n+i) = n0 + ampl*(cos(*(n+1)*(k0+k1)*double(i)) + 0.0*(-2.0*random_variable+1.0));
		//(exp(-4.0*log(10.0)*pow((double(i)-double(size)/2)/double(size),2.0)) - 0.1)*
	}*/

//	cout << cos(2.0*k0*n0*(size+4)) << endl;
	return;
}

int main()
{

	FILE * file_n = fopen("file_n.txt","w");
	FILE * file_E_plus = fopen("file_E_plus.txt","w");
	FILE * file_E_minus = fopen("file_E_minus.txt","w");
	FILE * file_Energy = fopen("file_Energy.txt","w");
	FILE * file_Efull = fopen("file_Efull.txt","w");
	FILE * file_diff = fopen("file_diff.txt","w");


	double * n = new double [size * sizeof(double)];
	set_n(n, file_n);
	double nr, nl;

	complex<double> EInit_plus (1.0, 0.0); 		//Волна, которая вылетает справа из среды 
	complex<double> EInit_minus (0.0, 0.0);		//Волна, которая влетает в среду справа (нулевая)

	complex<double> ELInit_plus (0.0, 0.0); 		//Волна, которая вылетает справа из среды 
	complex<double> ELInit_minus (1.0, 0.0);		//Волна, которая влетает в среду справа (нулевая)


	complex<double> * E_plus = new complex<double> [size+1]; 	// Массив амплитуд волн, которые распространяются слева направо
	*(E_plus+size) = EInit_plus;
	complex<double> * E_minus = new complex<double> [size+1];	// Массив амплитуд волн, которые распространяются справа налево
	*(E_minus+size) = EInit_minus;

	complex<double> * EL_plus = new complex<double> [size+1]; 	// Массив амплитуд волн, которые распространяются слева направо
	*(EL_plus) = ELInit_plus;
	complex<double> * EL_minus = new complex<double> [size+1];	// Массив амплитуд волн, которые распространяются справа налево
	*(EL_minus) = ELInit_minus;

	complex<double> * E_full = new complex<double> [size+1];	// Массив полного поля


	complex<double> Eleft_plus;
	complex<double> Eleft_minus;
	complex<double> Eright_plus;
	complex<double> Eright_minus;

	Eright_plus = EInit_plus;
	Eright_minus = EInit_minus;

	for (int i = size; i >= 1; i--)	//волна справа налево
	{
		nr = *(n+i);	// Текущие показатели преломления
		nl = *(n+i-1);

		complex<double> phaseR (0.0, k0*nr*double(i-size/2)); // i * kR * z
		complex<double> phaseL (0.0, k0*nl*double(i-size/2)); // 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 == size-1)
		//	cout << "E plus: " << Eright_plus  << "   Eright_minus: " << Eright_minus << "   exp(phaseR-L) " << exp(phaseR-phaseL) << "   exp(phaseR+L) " << exp(phaseL+phaseR) << endl;

		//Заполнение и сохранение массивов
		*(E_plus+i-1) = Eleft_plus;
		*(E_minus+i-1) = Eleft_minus;

		// Переход на следующую границу
		Eright_plus = Eleft_plus;
		Eright_minus = Eleft_minus;
	}

	Eleft_plus = ELInit_plus;
	Eleft_minus = ELInit_minus;

	for (int i = 1; i<= size; i++) //волна слева направо
	{
		nr = *(n+i);	// Текущие показатели преломления
		nl = *(n+i-1);

		complex<double> phaseR (0.0, k0*nr*double(i-size/2-1)); // i * kR * z
		complex<double> phaseL (0.0, k0*nl*double(i-size/2-1)); // 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 == 2)
		//	cout << "E plus: " << Eleft_plus  << "   E minus: " << Eleft_minus << "   exp(phaseR-L) " << exp(phaseR-phaseL) << "   exp(phaseR+L) " << exp(phaseL+phaseR) << endl;

		*(EL_plus+i) = Eright_plus;
		*(EL_minus+i) = Eright_minus;

		Eleft_plus = Eright_plus;
		Eleft_minus = Eright_minus;
	}

	for (int i = size; i >= 0; i--)
	{	
		*(E_minus+i) = *(E_minus+i) / *(E_plus);
		*(E_plus+i) = *(E_plus+i) / *(E_plus);
	}

		for (int i = 0; i <= size; i++)
	{
		*(EL_plus+i) = *(EL_plus+i) / *(EL_minus+size);
		*(EL_minus+i) = *(EL_minus+i) / *(EL_minus+size);
	}
	fprintf(file_E_plus,  "position    real E_plus   imag E_plus real EL_plus   imag EL_plus\n"); 
	for (int i = 0; i <= size; i++)
	{
			fprintf(file_n, "%e %e\n", double(i)/discretSize, *(n+i));
			fprintf(file_E_plus,  "%e %e %e %e %e\n", double(i)/discretSize,  real(*(E_plus+i)), imag(*(E_plus+i)),  real(*(EL_plus+i)), imag(*(EL_plus+i))  );
			fprintf(file_E_minus, "%e %e %e %e %e\n", double(i)/discretSize,  real(*(E_minus+i)),imag(*(E_minus+i)), real(*(EL_minus+i)), imag(*(EL_minus+i)) );
			//fprintf(file_diff, "%e %e %e\n", double(i)/discretSize, abs(*(E_plus+i)) - abs(*(EL_minus+size-i)), abs(*(EL_plus+i)) - abs(*(E_minus+size-i)));

			// file_Energy хранит вектор Пойнтинга	
			fprintf(file_Energy, "%e %e %e\n", double(i)/discretSize, ((pow(abs(*(E_plus+i)),2.0)) - pow(abs(*(E_minus+i)),2.0))*(*(n+i)), ((pow(abs(*(EL_plus+i)),2.0)) - pow(abs(*(EL_minus+i)),2.0))*(*(n+i)));

			//fprintf(file_Efull, "%e %e %e\n", double(i)/discretSize, arg( *(E_plus+i)), arg ( *(E_minus+size-i)) );
			fprintf(file_Efull, "%e %e %e\n", double(i)/discretSize, abs (*(E_plus+i) + (*(EL_plus+i)) ), abs( *(E_minus+i) + (*(EL_minus+i))) );
			//fprintf(file_Efull, "%e %e %e\n", double(i)/discretSize, abs (*(E_plus+i) + (*(E_minus+size-i)) ), abs( *(E_minus+i) + (*(E_plus+size-i))) );
	}

	fclose(file_n);
	fclose(file_E_plus);
	fclose(file_E_minus);
	fclose(file_Energy);
	fclose(file_Efull);
	fclose(file_diff);


	delete [] E_plus;
	delete [] E_minus;
	delete [] E_full;

	return 0;
}