One Hat Cyber Team

View File Name : find_period.cpp

#include <stdlib.h>

#include <iostream.h>

#include <fstream.h>

#include <string.h>

#include <math.h>

#include <stdio.h>







void Int2Str (int Num, char* Str, int Len)    

{    

  int NumTmp = Num, Tmp;    

  register int i;    

  strcpy(Str, "");    

  char TmpStr[100];    

    

  if(NumTmp < 0) return;    

    

  if(!NumTmp){    

    for(i = 0; i < Len; i++ )   strcat(Str,"0");    

    return;    

  }    

  while(NumTmp){    

    Tmp = NumTmp % 10;    

    switch(Tmp){    

    case 0: strcat(Str,"0"); break;    

    case 1: strcat(Str,"1"); break;    

    case 2: strcat(Str,"2"); break;    

    case 3: strcat(Str,"3"); break;    

    case 4: strcat(Str,"4"); break;    

    case 5: strcat(Str,"5"); break;    

    case 6: strcat(Str,"6"); break;    

    case 7: strcat(Str,"7"); break;    

    case 8: strcat(Str,"8"); break;    

    case 9: strcat(Str,"9"); break;    

    }    

    NumTmp /= 10;     

  }    

  NumTmp = strlen(Str);    

  for(i = 0; i < Len-NumTmp; i++ )   strcat(Str,"0");    

  NumTmp = strlen(Str);    

  strcpy(TmpStr, Str);    

  strcpy(Str, "");    

  for(i = NumTmp; i >= 0; i--) strncat(Str, TmpStr + i,1);    

  return;    

}    











main(int argn, char** argc)

{

  double x,y,z;

  double a=0.2, b=0.2, c=5.7, K=0.2;

  double TauFrom, TauTo, TauStep, Tau;

  double Time, TimePrev, x0=0.,  y0=-6.78,  z0=0.02;

  double *ArrXIni, *ArrYIni, *ArrZIni;

  double *ArrX, *ArrY, *ArrZ, *ArrX1, *ArrY1, *ArrZ1, *TmpAdr;

  unsigned long NumStepPer, NumStepPerMax, NumPerSkip, NumPerInt;

  register unsigned long i,j;

  unsigned long TauNumStep, NumPerWrite, NumPer;

  double  StepInt, Err = 1.0e-5, ErrStep, Diffr, StepWrite = 0.001;

  char Str[1000], fname[255], StrNum[10];

  FILE *FileOut, *FileTrj;



  if(argn < 11){

    cout << "Program to find  a limit cycle by its stabilization\n";

    cout << "Usage:  " << argc[0] << "  PerFrom  PerTo  PerStep  NumPerSkip  NumPerInt  StepIntgr ParK  File(out) NumPerPrn StepPrn\n";

    exit(1);

  }



  TauFrom = strtod(argc[1], (char**)NULL);

  TauTo = strtod(argc[2], (char**)NULL);

  if(TauFrom > TauTo) {cout << "Problem: TauFrom > TauTo. Exit.\n"; exit(1);}

  TauStep = strtod(argc[3], (char**)NULL);

  TauNumStep =  (unsigned long)((TauTo - TauFrom) / TauStep);

  NumPerSkip = atoi(argc[4]);

  NumPerInt = atoi(argc[5]);

  StepInt = strtod(argc[6], (char**)NULL);

  ErrStep = StepInt / 2.;

  K = strtod(argc[7], (char**)NULL);

  strcpy(fname, argc[8]);

  NumPerWrite = atoi(argc[9]);

  StepWrite = strtod(argc[10], (char**)NULL);



  strcpy(Str, fname);

  strcat(Str,".diz");

  ofstream ToDiz(Str);

  if(!ToDiz) {

    cout << "Cannot open file " << Str << ". Exit. \n" << flush;

    exit(1);

  }

  

  ToDiz << "Program to find  a limit cycle by its stabilization\n" << flush;

  ToDiz << "Range of periods to be checked:\n";

  ToDiz << "From = " << TauFrom << ",  To = " << TauTo << ",  Step = " <<  TauStep << 

                                      ", NumSteps = " << TauNumStep +1  << endl << flush;

  ToDiz << "Integration will be done over  " << NumPerInt <<  " periods" << endl;

  ToDiz << "To Skip : " << NumPerSkip << " periods.\n";

  ToDiz << "Step of integration: " << StepInt << endl << flush;

  ToDiz << "Parameters of the System: a = " << a << ", b = " << b << ", c = " << c << endl << flush;

  ToDiz << "Force of Delay = " << K << endl << endl;

  ToDiz << "The last  " << NumPerWrite << " period will be saved to " << fname <<".txxxxx\n" ;   

  ToDiz << "Step To Write:  " << StepWrite << endl << flush;

    

  cout << "Range of periods to be checked:\n";

  cout << "From = " << TauFrom << ",  To = " << TauTo << ",  Step = " <<  TauStep << 

                                      ", NumSteps = " << TauNumStep +1  << endl << flush;

  cout << "Integration will be done over  " << NumPerInt <<  " periods" << endl;

  cout << "To Skip : " << NumPerSkip << " periods.\n";

  cout << "Step of integration: " << StepInt << endl << flush;

  cout << "Parameters of the System: a = " << a << ", b = " << b << ", c = " << c << endl << flush;

  cout << "Force of Delay = " << K << endl;

  cout << "The last  " << NumPerWrite << " period will be saved to " << fname <<".txxxxx\n" ;   

  cout << "Step To Write:  " << StepWrite << endl << flush;

  

  ToDiz << "Doing Transient of " << NumPerSkip << "  max. periods = " << 

                              NumPerSkip*TauTo << " time... -> " << flush;

/*  cout << "Doing Transient of " << NumPerSkip << "  max. periods = " << 

                              NumPerSkip*TauTo << " time... -> " << flush; */

  NumStepPerMax = (unsigned long)(TauTo / StepInt);

   if(NumStepPerMax*StepInt < TauTo) NumStepPerMax++;



  ArrXIni =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  ArrYIni =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  ArrZIni =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  if(!ArrXIni || !ArrYIni || !ArrZIni) {

    cout << "Cannot allocate memory for StateIni arrays. Exit. \n";   exit(1);

  }



  for(i = 0; i <= NumPerSkip; i++) {

/*    cout << "Period:  " << i << endl; */

    if( i == NumPerSkip) {ArrXIni[0] = x0; ArrYIni[0] = y0; ArrZIni[0] = z0;}

    for(j = 1; j <= NumStepPerMax; j++) {

      if( i == NumPerSkip) {ArrXIni[j] = x0; ArrYIni[j] = y0; ArrZIni[j] = z0;}

      x=x0+StepInt*(-y0-z0);

      y=y0+StepInt*(x0+a*y0);

      z=z0+StepInt*(b+z0*(x0-c));

      x0=x;  y0=y;  z0=z;

    }    

  }  

  ArrXIni[NumStepPerMax+1] = x; 

  ArrYIni[NumStepPerMax+1] = y; 

  ArrZIni[NumStepPerMax+1] = z; 

  

  ToDiz << "Done! \n" << flush;

/*  cout << "Done! \n" << flush; */



  ArrX =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  ArrY =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  ArrZ =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  ArrX1 =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  ArrY1 =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  ArrZ1 =(double*)malloc((NumStepPerMax+2)*sizeof(double));

  if(!ArrX1 || !ArrY1 || !ArrZ1 ||  !ArrX1 || !ArrY1 || !ArrZ1) {

    cout << "Cannot allocate memory for State arrays. Exit. \n";   exit(1);

  }

  

  ToDiz << "\nStart to look for Limit Cycles \n"  << flush; 



  FileOut = fopen(fname, "wt");

  

  if(!FileOut) { ToDiz << "Cannot open file. " << fname << " Exit." << endl;    exit(1); }

  

  NumPer = 0;

  for(Tau = TauFrom; Tau <= TauTo; Tau += TauStep){

    NumPer++;

    NumStepPer = (unsigned long)(Tau / StepInt);        

    if(NumStepPer*StepInt < Tau) NumStepPer++;

    ToDiz << "N " << NumPer << " :  Period to test:  " << Tau << ",  NumSteps =  " << NumStepPer; 

    if(NumPerWrite) {

      strcpy(Str, fname);

      Int2Str(NumPer, StrNum,5);

      strcat(Str,".t");   strcat(Str, StrNum);

      FileTrj = fopen(Str,"wt");

      if(!FileTrj) {ToDiz << "Cannot open file " << Str << ". Exit. " << endl << flush; exit(1);}

      else ToDiz << "  File:  " << StrNum  << endl << flush;  

    }

    else  ToDiz  << endl << flush;  



    TimePrev = Time = 0;



    for(i = 0;  i <= NumStepPer; i++)

      ArrX[i] = ArrXIni[i]; ArrY[i] = ArrYIni[i]; ArrZ[i] = ArrZIni[i];

    Diffr = 0;

    for(i = 0; i <= NumPerInt; i++){

      ToDiz << "               Skip period:  " << i  << endl << flush;

      ArrX1[0] = ArrX[NumStepPer];  

      ArrY1[0] = ArrY[NumStepPer];  

      ArrZ1[0] = ArrZ[NumStepPer];



      if(i >  NumPerInt - NumPerWrite) 

        fprintf(FileTrj,"%15.8lf  %15.8lf  %15.8lf %15.8lf\n", Time, ArrX1[0], ArrY1[0], ArrZ1[0]);



      for(j = 1; j <= NumStepPer; j++) {



        x0 = ArrX1[j-1]; y0 = ArrY1[j-1]; z0 = ArrZ1[j-1];

        ArrX1[j]= x0+StepInt*((-y0-z0)+K*(ArrX[j-1]-x0));

        ArrY1[j]=y0+StepInt*((x0+a*y0)+K*(ArrY[j-1] - y0));

        ArrZ1[j]=z0+StepInt*((b+z0*(x0-c))+K*(ArrZ[j-1]-z0));

        if( i == NumPerInt) Diffr += fabs(ArrX[j-1] - x0)+fabs(ArrY[j-1] - y0)+fabs(ArrZ[j-1] - z0);

        Time = i*Tau + j*StepInt;

        if((i >  NumPerInt - NumPerWrite)) { 

          fprintf(FileTrj,"%15.8lf  %15.8lf  %15.8lf %15.8lf\n ", Time, ArrX1[j], ArrY1[j], ArrZ1[j]);

          fflush(FileTrj);

          TimePrev = Time; 

        }

      }

      TmpAdr = ArrX; ArrX = ArrX1; ArrX1 = TmpAdr;

      TmpAdr = ArrY; ArrY = ArrY1; ArrY1 = TmpAdr;

      TmpAdr = ArrZ; ArrZ = ArrZ1; ArrZ1 = TmpAdr;

      if(i >  NumPerInt - NumPerWrite && FileTrj) fprintf(FileTrj, "\n");

    }

    if(NumPerWrite && FileTrj) fclose(FileTrj);

    fprintf(FileOut,"%12.8lf   %12.8lf\n", Tau, log(Diffr *= StepInt));

    fflush(FileOut);

  }





  ToDiz << "Calculations are finished\n";

  fclose( FileOut );

  free(ArrXIni);

  free(ArrYIni);

  free(ArrZIni);

  free(ArrX);

  free(ArrY);

  free(ArrZ);

  free(ArrX1);

  free(ArrY1);

  free(ArrZ1);



}