Puncare_Map_ADgyro_v.1.1_win.cpp

//-------------------------------------------------------------------------------------------
//              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 <math.h>
#include "time_measurement.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 M_PIl //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

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


//------------------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, double y[DIM])//Puancare cross-section formula
{
	double value;
	value=sin((nu*time)/2);// repetition of the phase:  nu*time=nu*time+2*Pi
	return value;
}



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

//----------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
{
//Reference of extern global variables:

	//local variebles:
	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_name(FILE *f_ini, double *x, char *name_file){
	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",name_file);
}

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

void process_trajectory(FILE *f_ini, FILE* f_lL_res, int k, int NumPointsOnMap, double *_StepIntegr, double G, double TolIntegr, double TolOfSection){
	double StepIntegr = *_StepIntegr;
	double t = 0;
	double x[DIM], X[DIM];//arrays for start and final point of one iteration
	double T;
	int NumPointsOnMap_NOW = 0;
	char name_file[100];
	double rk_temp, StepTemp;
	int i;
	FILE *f_res;

	printf("%4d",k );
	read_x_name(f_ini, x, name_file);

	f_res=fopen(name_file,"w");
	//end of read of line with initial values for k trajectory

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

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

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

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

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



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

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

				StepIntegr=StepTemp;//restor step size

				fprintf(f_lL_res,"%15.10f %15.10f\n",fmod((fmod(x[0],2*PI)+2*PI),2*PI),x[1]/G);//saving of results to the common file lL.txt
			}// 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

	fclose(f_res);
	*_StepIntegr = StepIntegr;
}

//  MAIN PROGRAM BODY
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);

	for(k=1;k<=NumTrajectories;k++)
		process_trajectory(f_ini, f_lL_res, k, NumPointsOnMap, &StepIntegr, G, TolIntegr, TolOfSection);

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

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