Puncare_Map_ADgyro_v.1.1_win_parallel.c

//-------------------------------------------------------------------------------------------
//              This file is the MAIN part of the program complex "Poincare_Map"
//              Author: Anton Doroshin <doran@inbox.ru>
//              Samara State Aerospace University (Russia)
//              Version 1.1 (04.04.2008)
//				Version 2.0 (07.05.2013)
//
//              THIS IS FREE SOFTWARE on GNU GPL BASE
//
//              YOU MAY USE, MODIFY, DISTRIBUTE THIS PROGRAM
//
//--------------------------------------------------------------------------------------------
//      Program for Poincare map of dinamical systems evaluating (On the Gyrostat's example base).
//      Win-Lin multiplatform program
//      Program have been tested on UBUNTU and WINDOWS_XP platform
//
//		WIN_DEVELOP TOOL:
//      Microsoft Visual Studio 2008  Version 9.0.21022.8 RTM
//      Microsoft .NET Framework Version 3.5
//      Installed Edition: VC Express
//      Microsoft Visual C++ 2008   91909-152-0000052-60150
//--------------------------------------------------------------------------------------------

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "time_measurement.h"
#include <omp.h>

#define DIM 2 //Dimension of the dynamical system <--second order as an example
#define EGOROVS_CONST 8 //Egolov's constant for integration step control
#define PI 3.1415926535897932384626433832795029L //Pi-magic-number is the base of the Universe

//some constants take place:
#define A1		5.0
#define A2		15.0
#define B2		8.0
#define C1		4.0
#define C2		7.0
#define nu		1.0
#define eps		0.60
#define DELTA	3.0

struct trajectory_processing_input_data{
	double x[DIM];
	double NumPointsOnMap;
	double StepIntegr;
	double TolIntegr;
	double TolOfSection;
	double G;
};

struct trajectory_processing_output_data{
	double *points;
};

//------------------------------------------------------------------------------------------


//------------------RIGHT-HANDS OF DIFFERENTIAL EQUATIONS SYSTEM----------------------------
//
//    Following notations take place:
//    time - time or other independent variable
//    y[0] ... y[DIM-1] - phase coordinates of dynamical system
//
//						y[0] is the first coordinate
//
//------------------------------------------------------------------------------------------

double F(int number_of_equation, double time, double y[DIM], double G)//Right-hand of dynamical system
{
	double value;

	switch (number_of_equation)
	{ //   RIGHT-HANDS OF DIFFERENTIAL EQUATIONS SYSTEM:

/* If You have another dimention (DIM is define in "#define DIM YOU_DIM_NUMBER") of system then add or remove your equations*/

	case 0: /* First equation */
		value = y[1]*(1/C2-cos(y[0])*cos(y[0])/(A1+B2)-sin(y[0])*sin(y[0])/(A1+A2))-DELTA/C2-eps/C2/nu*sin(nu*time);
		break;
	case 1: /* 2_equation */
		value = -(G*G-y[1]*y[1])*(1/(A1+A2)-1/(A1+B2))*sin(y[0])*cos(y[0]);
		break; //
	default:
		value=0;
		break;
	}

	return value;
}// 


//--------------------POINCARE SURFACE (CROSSECTION CONDITION) -----------------------------------------
double Section(double time)//Puancare cross-section formula
{
	return sin((nu*time)/2);// repetition of the phase:  nu*time=nu*time+2*Pi
}



//-----------------------------------------------------------------------------

//----------INTEGRATION METHOD-------------------------------------------------
//
//Function RK45() return the magnitude of Egorov's control term
//
//
//You also can change the method-formulas
//
//Following notations take place:
//DIM - dimention of dymamical system
//t - time
//X[0] ... X[DIM-1] - phase coordinates of dynamical system
//k1[i] ... k4[i] - arrays of coefficients Classical_RK45-Method
//-----------------------------------------------------------------------------

double RK45(double t, double *T, double *x, double *X, double StepIntegr, double G)//Function for evaluation of arrays kj[], X[]
			//Function return control term value
{
	double k1[DIM],k2[DIM],k3[DIM],k4[DIM]; //arrays for RK45 and ControlTerm evaluation
	double x_k1[DIM],x_k2[DIM],x_k3[DIM];
	int i;
	double value=0;
	double temp=0;

	for (i=0; i<DIM; i++)
	{
		k1[i]=StepIntegr*F(i, t, x, G);
		x_k1[i]=x[i]+0.5*k1[i];
	}

	for (i=0; i<DIM; i++)
	{
		k2[i]=StepIntegr*F(i,t+0.5*StepIntegr,x_k1, G);
		x_k2[i]=x[i]+0.5*k2[i];
	}

	for (i=0; i<DIM; i++)
	{
		k3[i]=StepIntegr*F(i,t+0.5*StepIntegr,x_k2, G);
		x_k3[i]=x[i]+k3[i];
	}

	for (i=0; i<DIM; i++)
	{
		k4[i]=StepIntegr*F(i,t+StepIntegr,x_k3, G);
	}
	
	for (i=0; i<DIM; i++)
	{
		X[i]=x[i]+(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6;
	}
	*T=t+StepIntegr;

	for (i=0;i<DIM;i++)//evaluation of Egorov's control term:
	{
		temp=fabs(k1[i]-k2[i]-k3[i]+k4[i])*2/3;
		if (temp>value) {
			value=temp;
		}
	}
	return value;
}
//----------------------------------------------------------------------------

void read_x(FILE *f_ini, double *x){
	char temp[500];
	int i;

	fscanf(f_ini,"%s",temp);
	fscanf(f_ini,"%s",temp);

	int bufint;
	fscanf(f_ini,"%d",&bufint);
	for (i=0;i<DIM;i++)
	{
		fscanf(f_ini,"%lf", &x[i]);
	}
	fscanf(f_ini,"%s",temp);
}

void read_some_params(FILE *f_ini, int *NumTrajectories, int *NumPointsOnMap, double *StepIntegr, double *TolIntegr, double *TolOfSection, double *G){
	char temp[500];
	fscanf(f_ini,"%s",temp);

	fscanf(f_ini,"%s",temp);

	fscanf(f_ini,"%s",temp);

	fscanf(f_ini,"%d %d %lf %lf %lf\n",NumTrajectories, NumPointsOnMap, StepIntegr, TolIntegr, TolOfSection);

	fscanf(f_ini,"%s",temp);
	fscanf(f_ini,"%s",temp);
	fscanf(f_ini,"%lf  \n", G);
}

struct trajectory_processing_output_data process_trajectory(struct trajectory_processing_input_data input){
	double t = 0;
	double x[DIM], X[DIM];//arrays for start and final point of one iteration
	double T;
	int NumPointsOnMap_NOW = 0;
	double rk_temp, StepTemp;
	double StepIntegr = input.StepIntegr;
	int i;
	struct trajectory_processing_output_data res;
	res.points = (double*)malloc(input.NumPointsOnMap*sizeof(*res.points)*DIM);

	for(i=0;i<DIM;i++)
		x[i]=input.x[i];

	while (NumPointsOnMap_NOW<input.NumPointsOnMap) //finding of each point of the Poincare-map
	{//1
		//---------------ONE ITERATION WITH CONTROL TERM---------------

		rk_temp=RK45(t, &T, x, X, StepIntegr, input.G);//execute of one iteration of Runge-Kutta method

		while (rk_temp>input.TolIntegr)
		{//autochanging of StepIntegr_SIZE
			StepIntegr/=2;
			rk_temp=RK45(t, &T, x, X, StepIntegr, input.G);
		}

		while (rk_temp<=input.TolIntegr/EGOROVS_CONST)
		{
			StepIntegr*=2;
			rk_temp=RK45(t, &T, x, X, StepIntegr, input.G);
		}//autochanging of StepIntegr_SIZE */

		//----------------END OF ITERATION-----------------------------



		//-----------INTERSECTION FINDING------------------------------
		double SectionBefore, SectionAfter;// values of section function befor and after intersection
		SectionBefore=Section(t);
		SectionAfter=Section(T);

		if (SectionBefore*SectionAfter<0)//then we have intersection on this iteration-interval
		{//2
			if (SectionBefore>0)//(X[3]>x[3])
			{// direction of trajectory transpotation control
				StepTemp=StepIntegr;//Let's save in memory step size befor finding

				while (fabs(SectionAfter)>input.TolOfSection) //correction of tolerance of point
				{//3
					if (SectionBefore*SectionAfter<0) StepIntegr*=0.5; else StepIntegr*=1.5;
					RK45(t, &T, x, X, StepIntegr, input.G);
					SectionAfter=Section(T);
					//This method corespond to "half-secant" ideology
				}//3

				StepIntegr=StepTemp;//restor step size

				res.points[NumPointsOnMap_NOW*DIM] = fmod((fmod(x[0],2*PI)+2*PI),2*PI);
				res.points[NumPointsOnMap_NOW*DIM+1] = x[1]/input.G;
				NumPointsOnMap_NOW++;

			}// direction of trajectory transpotation control
		}//2
		//----CHANGING OF START-FINAL POINTS FOR FURTHER ITERATIONS-----
		for (i=0;i<DIM;i++)
			x[i]=X[i];
		t=T;
	}//1
	return res;
}

int main(int argc, char* argv[])
{
	FILE *f_ini, *f_lL_res;
	int NumTrajectories;
	double StepIntegr;
	double TolIntegr, TolOfSection;
	int NumPointsOnMap; //number of points on Puancare map for one trajectory
	double G;
	int k;
	unsigned long int start_time;

	f_lL_res	= fopen("./RESULT/lL.txt","w");
	f_ini		= fopen("ini.txt","r");

	printf("I'm program complex \"Puancare_Map\".\n" );
	printf("I'm working, please wait.\n" );

	start_time = get_current_time_in_us();

	read_some_params(f_ini, &NumTrajectories, &NumPointsOnMap, &StepIntegr, &TolIntegr, &TolOfSection, &G);

#pragma omp parallel for schedule(dynamic)
	for(k=1;k<=NumTrajectories;k++){
		struct trajectory_processing_input_data input;
		struct trajectory_processing_output_data output;
		printf("Processing trajectrory %d\n",k );
		input.NumPointsOnMap = NumPointsOnMap;
		input.StepIntegr = StepIntegr;
		input.G = G;
		input.TolIntegr = TolIntegr;
		input.TolOfSection = TolOfSection;

#pragma omp critical(datainput)
		read_x(f_ini, input.x);

		output = process_trajectory(input);

#pragma omp critical(dataoutput)
		{
			for(int i=0; i<NumPointsOnMap; i++){
				fprintf(f_lL_res,"%15.10f %15.10f\n",output.points[i*DIM],output.points[i*DIM+1]);
			}
		}
		free(output.points);
	}

	printf("\nDone in: %ld\n", get_current_time_in_us() - start_time);

	fclose(f_ini);
	fclose(f_lL_res);
	return 0;
}