refract.cc

#include <iostream>
#include <cmath>
#include <stdio.h>
using namespace std;


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");


	int size = 100000;
	double ampl = 5.0e-1;

	double E_init_plus_Re = 1.0;
	double E_init_plus_Im = 0.0;
	double E_init_minus_Re = 0.0;
	double E_init_minus_Im = 0.0;

	double E_left_plus_Re; //             ___________________________\  |
	double E_left_plus_Im; //                                        /  |

	double E_left_minus_Re; //           /____________________________  |
	double E_left_minus_Im; //           \                              |    

	//double E_right_plus_Re; //           | _____________________________\ 
	//double E_right_plus_Im; //           |                              / 

	//double E_right_minus_Re; //          | /_____________________________ 
	//double E_right_minus_Im; //          | \                               

	double * n = new double [size * sizeof(double)];
	double * E_plus_Re = new double [size * sizeof(double)];
	double * E_plus_Im = new double [size * sizeof(double)];
	double * E_minus_Re = new double [size * sizeof(double)];
	double * E_minus_Im = new double [size * sizeof(double)];

	for(int i = 0; i < size; i++){
		*(n+i) = 1.0 + ampl*(cos(1/100.0*2*2*3.1415926/100*i)+1.0);
		/*if (i > size/2)
			*(n+i) = 2.0;
		else
			*(n+i) = 1.0;*/
		/*if ((i/100)%2 == 0)
		{	
			*(n+i) = sqrt(2.0);
		}
		else	
			*(n+i) = 1.0;*/
		fprintf(file_n, "%e\n", *(n+i));
		//cout << i << endl;
	}

	double E_right_plus_Re = E_init_plus_Re;
	double E_right_plus_Im = E_init_plus_Im;
	double E_right_minus_Re = E_init_minus_Re;
	double E_right_minus_Im = E_init_minus_Im;

	double P0, P_plus, P_minus;

	for (int i = size-1; i >= 0; i--)
	{
		//if( i < size/100.0 || i > size*0.99){
		//if ((i+1.0)/100.0 == int((i+1)/100.0)){
			//fprintf(file_E_plus, "%e\n", sqrt(pow(E_right_plus_Re,2.0) + pow(E_right_plus_Im,2.0)));
		//	fprintf(file_E_plus, "%e %e %e\n", E_right_plus_Re, E_right_plus_Im, sqrt(pow(E_right_plus_Re,2.0) + pow(E_right_plus_Im,2.0)));
		//	fprintf(file_E_minus, "%e %e %e\n", E_right_minus_Re, E_right_minus_Im, sqrt(pow(E_right_minus_Re,2.0) + pow(E_right_minus_Im,2.0)));
		//}

		*(E_plus_Re+i) = E_right_plus_Re;
		*(E_plus_Im+i) = E_right_plus_Im;
		*(E_minus_Re+i) = E_right_minus_Re;
		*(E_minus_Im+i) = E_right_minus_Im;

		double nr = *(n+i);
		double nl;
		if (i == 0)
			nl = nr;
		else
			nl = *(n+i-1);

		//cout << nl << "  " << nr << endl;
		//cin.get();
	
		E_left_plus_Re = (cos(2*3.1415926*(-nl/100 + (nr-nl)*(i+1)/100.0))*E_right_plus_Re  - sin(2*3.1415926*(-nl/100 + (nr-nl)*(i+1)/100.0))*E_right_plus_Im)*(nl+nr)/(2*nl);
		E_left_plus_Re +=(cos(2*3.1415926*(-nl/100 - (nr+nl)*(i+1)/100.0))*E_right_minus_Re - sin(2*3.1415926*(-nl/100 - (nr+nl)*(i+1)/100.0))*E_right_minus_Im)*(nl-nr)/(2*nl);

		E_left_plus_Im = (sin(2*3.1415926*(-nl/100 + (nr-nl)*(i+1)/100.0))*E_right_plus_Re  + cos(2*3.1415926*(-nl/100 + (nr-nl)*(i+1)/100.0))*E_right_plus_Im)*(nl+nr)/(2*nl);
		E_left_plus_Im +=(sin(2*3.1415926*(-nl/100 - (nr+nl)*(i+1)/100.0))*E_right_minus_Re + cos(2*3.1415926*(-nl/100 - (nr+nl)*(i+1)/100.0))*E_right_minus_Im)*(nl-nr)/(2*nl);

		E_left_minus_Re = (cos(2*3.1415926*(nl/100 + (nl+nr)*(i+1)/100.0))*E_right_plus_Re  - sin(2*3.1415926*(nl/100 + (nl+nr)*(i+1)/100.0))*E_right_plus_Im)*(nl-nr)/(2*nl);
		E_left_minus_Re +=(cos(2*3.1415926*(nl/100 + (nl-nr)*(i+1)/100.0))*E_right_minus_Re - sin(2*3.1415926*(nl/100 + (nl-nr)*(i+1)/100.0))*E_right_minus_Im)*(nl+nr)/(2*nl);

		E_left_minus_Im = (sin(2*3.1415926*(nl/100 + (nl+nr)*(i+1)/100.0))*E_right_plus_Re  + cos(2*3.1415926*(nl/100 + (nl+nr)*(i+1)/100.0))*E_right_plus_Im)*(nl-nr)/(2*nl);
		E_left_minus_Im +=(sin(2*3.1415926*(nl/100 + (nl-nr)*(i+1)/100.0))*E_right_minus_Re + cos(2*3.1415926*(nl/100 + (nl-nr)*(i+1)/100.0))*E_right_minus_Im)*(nl+nr)/(2*nl);

		E_right_plus_Re = E_left_plus_Re;
		E_right_plus_Im = E_left_plus_Im;
		E_right_minus_Re = E_left_minus_Re;
		E_right_minus_Im = E_left_minus_Im;
	}

	double ampl_init = sqrt(pow(*(E_plus_Re),2.0) + pow(*(E_plus_Im),2.0));
	double sin_init = *(E_plus_Im)/ampl_init;
	double cos_init = *(E_plus_Re)/ampl_init;

	for (int i = 0; i < size; i++)
	{
		E_right_plus_Re = *(E_plus_Re+i);
		E_right_plus_Im = *(E_plus_Im+i);
		E_right_minus_Re = *(E_minus_Re+i);
		E_right_minus_Im = *(E_minus_Im+i);

		*(E_plus_Re+i) = (E_right_plus_Re*cos_init + E_right_plus_Im*sin_init)/ampl_init;
		*(E_plus_Im+i) = (E_right_plus_Re*(-sin_init) + E_right_plus_Im*cos_init)/ampl_init;
		*(E_minus_Re+i) = (E_right_minus_Re*cos_init + E_right_minus_Im*sin_init)/ampl_init;
		*(E_minus_Im+i) = (E_right_minus_Re*(-sin_init) + E_right_minus_Im*cos_init)/ampl_init;
		
		//if (i < 0.01*size || i > 0.99*size){
		/*	
			*/
		//}
			fprintf(file_E_plus, "%e %e %e\n", *(E_plus_Re+i), *(E_plus_Im+i), sqrt(pow(*(E_plus_Re+i),2.0) + pow(*(E_plus_Im+i),2.0)));
			fprintf(file_E_minus, "%e %e %e\n", *(E_minus_Re+i), *(E_minus_Im+i), sqrt(pow(*(E_minus_Re+i),2.0) + pow(*(E_minus_Im+i),2.0)));

			P_plus = (pow(*(E_plus_Re+i),2.0) + pow(*(E_plus_Im+i),2.0)) * (*(n+i));
			P_minus = (pow(*(E_minus_Re+i),2.0) + pow(*(E_minus_Im+i),2.0)) * (*(n+i));
			P0 = P_plus - P_minus;
			fprintf(file_Energy, "%e\n", P0);
}

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

	return 0;
}