Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
//-------------------------------------------------------------------------------------------
// 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"
#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);
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;
read_x(f_ini, input.x);
output = process_trajectory(input);
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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.