#include <cstdlib>
#include <iostream>

#include <vector>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <conio.h>


#define MAX 5000000

using namespace std;


float B[MAX][4];

float y[MAX];


float C[100];
float AverageMean[100];
float MaxPeak[100];
float MinValley[100];

float MaxAmp[100];
float AverageAmp[100];

float ymax=0.;

float slope1,slope2;

float mina,maxa;

float t;
float scale;
float X,Y;

long ijk;

long kv;

long hold;
long i,j,k,n;
long num;

long nkv=0;

int ic,iscale;

long N,NP;

long last_a;

FILE *pFile[4];
char filename[4][FILENAME_MAX];

int main()
{

	printf("\n ");
	printf("\n rainflow.cpp  ver 2.0  November 3, 2011 \n");
	printf("\n by Tom Irvine  Email: tomirvine@aol.com \n");
	printf("\n ASTM E 1049-85 (2005) Rainflow Counting Method \n");

	printf("\n The input file must be a time history. \n");
	printf("\n Select format: ");
	printf("\n   1=amplitude ");
	printf("\n   2=time & amplitude \n");

	scanf("%d",&ic);

	if(ic==1)
	{
   		printf( "\n\n The base input file must contain one column: unit \n");
	}
	else
	{
   		printf( "\n\n The base input file must contain two columns: \n");
		printf( " time & unit    \n");
	}

	printf( "\n Input filename \n");
	scanf("%s",filename[0]);

    pFile[0] = fopen(filename[0], "rb");

	while(pFile[0] == NULL )
	{

		printf("\n Failed to open file: %s \n", filename[0]);

		printf("\n Please enter the input filename: \n");

		scanf("%s",filename[0]);

		pFile[0] = fopen(filename[0], "rb");

	}
	printf("\n File: %s opened. \n", filename[0]);

	printf( "\n\n Enter the output filename: \n");
	scanf("%s",filename[1]);

	pFile[1]=fopen(filename[1], "w");

    strcpy(filename[2],"rainflow_graph.out");
	pFile[2]=fopen(filename[2], "w");

	i=0;

	if(ic==1)
	{
		while( fscanf(pFile[0],"%f",&y[i])>0)
		{
			i++;

			if(i==MAX)
			{
				break;
			}
		}
	}
	else
	{
		while( fscanf(pFile[0],"%f %f",&t,&y[i])>0)
		{
			i++;

			if(i==MAX)
			{
				break;
			}
		}
	}
	NP=i+1;


//	printf("\n ref 1: last_a = %ld \n",last_a);

	printf("\n ");
	printf("\n Multiply data by scale factor?");
	printf("\n 1=yes  2=no \n");
	scanf("%d",&iscale);

	if(iscale==1)
	{
    	printf("\n Enter scale factor \n");
    	scanf("%f",&scale);

		for(i=0;i<NP;i++)
		{
			y[i]*=scale;
		}
	}


	 vector<float> a;


	k=0;
//	a[k]=y[k];
    a.push_back(y[k]);

	k=1;


    for(i=1;i<(NP-1);i++)
	{
    	if( y[i]>y[i-1] && y[i]>y[i+1])
		{
       // 	a[k]=y[i];
        	    a.push_back(y[i]);
        	k++;
		}
		else
		{
            if( y[i]<y[i-1] && y[i]<y[i+1])
            {
 //               a[k]=y[i];
         	    a.push_back(y[i]);
                k++;
            }
		}

	}

	last_a=k-1;

	hold=last_a;

//	a[k]=y[NP-1];
	a.push_back(y[NP-1]);

	mina=100000;
	maxa=-mina;

	for(i=0;i<=last_a;i++)
	{
		if(a[i]<mina)
		{
			mina=a[i];
		}
		if(a[i]>maxa)
		{
			maxa=a[i];
		}
	}

	num=long(maxa-mina)+1;

	n=0;
	i=0;
	j=1;

	float sum=0;

	kv=0;

    printf("\n percent completed \n");

	while(1)
	{
    	Y=(fabs(a[i]-a[i+1]));
    	X=(fabs(a[j]-a[j+1]));

    	if(X>=Y && Y>0)
		{
            if(Y>ymax)
            {ymax=Y;}

        	if(i==0)
			{

           		n=0;
           		sum+=0.5;
                B[kv][3]=a[i+1];
                B[kv][2]=a[i];
                B[kv][1]=0.5;
                B[kv][0]=Y;
                kv++;
/*
				for(k=0;k<last_a;k++)
				{
					a[k]=a[k+1];
				}
*/
				a.erase (a.begin());


				last_a--;

           		i=0;
           		j=1;
			}
        	else
			{
           		sum+=1;
                B[kv][3]=a[i+1];
                B[kv][2]=a[i];
                B[kv][1]=1.;
                B[kv][0]=Y;
                kv++;

 /*
            		for(k=0;k<i;k++)
				{
                        a.push_back(y[k]);
                    	n++;
				}
  		        for(k=i+2;k<=last_a;k++)
				{

                        a.push_back(y[k]);
                    	n++;
				}
	*/

           		n=0;

                a.erase (a.begin()+(i+1));
                a.erase (a.begin()+i);


				last_a-=2;

           		i=0;
           		j=1;
			}

			nkv++;

			if(nkv==3000)
			{
			   printf(" %8.4g \n",(sum/(hold/2))*100.);
			   nkv=0;
           }
		}
    	else
		{
        	i++;
        	j++;
		}

    	if((j+1)>last_a)
		{
        	break;
		}
	}


	for(i=0;i<(last_a);i++)
	{
    	Y=(fabs(a[i]-a[i+1]));
    	if(Y>0)
		{
        	sum+=0.5;
            B[kv][3]=a[i+1];
            B[kv][2]=a[i];
            B[kv][1]=0.5;
            B[kv][0]=Y;
            kv++;
		}
        if(Y>ymax)
        {ymax=Y;}
	}

double L[20];

L[1]=0;
L[2]=2.5;
L[3]=5;
L[4]=10;
L[5]=15;
L[6]=20;
L[7]=30;
L[8]=40;
L[9]=50;
L[10]=60;
L[11]=70;
L[12]=80;
L[13]=90;
L[14]=100;


for(ijk=1;ijk<=14;ijk++)
{
   L[ijk]*=ymax/100.;

   C[ijk]=0.;
   AverageMean[ijk]=0.;
   MaxPeak[ijk]=-1.0e+20;
   MinValley[ijk]= 1.0e+20;

   MaxAmp[ijk]=-1.0e+20;
   AverageAmp[ijk]= 1.0e+20;

//   printf("  L[%ld]=%g \n",ijk,L[ijk]);
}

kv--;

for(ijk=13;ijk>=0;ijk--)
{
   AverageAmp[ijk]=0.;
}


for(i=0;i<=kv;i++)
{
     Y=B[i][0];
//     printf("i=%d Y=%g \n",i,Y);

     for(ijk=13;ijk>=0;ijk--)
     {

        if(Y>=L[ijk] && Y<=L[ijk+1])
        {
            C[ijk]=C[ijk]+B[i][1];
            AverageMean[ijk]+=B[i][1]*(B[i][3]+B[i][2])*0.5; // weighted average

            if(B[i][3]>MaxPeak[ijk])
            {
               MaxPeak[ijk]=B[i][3];
            }
            if(B[i][2]>MaxPeak[ijk])
            {
               MaxPeak[ijk]=B[i][2];
            }

            if(B[i][3]<MinValley[ijk])
            {
               MinValley[ijk]=B[i][3];
            }
            if(B[i][2]<MinValley[ijk])
            {
               MinValley[ijk]=B[i][2];
            }

            if(Y>MaxAmp[ijk])
            {
               MaxAmp[ijk]=Y;
            }

            AverageAmp[ijk]+=B[i][1]*Y*0.5;

            break;
        }
     }
}

for(ijk=1;ijk<=14;ijk++)
{
   if(C[ijk]>0)
   {
      AverageMean[ijk]/=C[ijk];
       AverageAmp[ijk]/=C[ijk];
   }
   MaxAmp[ijk]/=2.;

   if(C[ijk]<0.5)
   {
       AverageAmp[ijk]=0.;
       MaxAmp[ijk]=0.;
       AverageMean[ijk]=0.;
       MinValley[ijk]=0.;
       MaxPeak[ijk]=0.;
   }
}

   printf("\n Amplitude = (peak-valley)/2 \n");

   fprintf(pFile[1],"\n Amplitude = (peak-valley)/2 \n");

//*****************************************************************************


	          printf("\n ");
	          printf("\n          Range            Cycle       Ave      Max     Ave     Min       Max");
	          printf("\n         (units)           Counts      Amp      Amp     Mean    Valley    Peak \n");

	for(i=13;i>=1;i--)
    {
    	printf("  %8.2lf to %8.2lf\t%8.1lf\t%6.4g\t%6.4g\t%6.4g\t%6.4g\t %6.4g\n",L[i],L[i+1],C[i],AverageAmp[i],MaxAmp[i],AverageMean[i],MinValley[i],MaxPeak[i]);
    }

	fprintf(pFile[1],"\n ");
	fprintf(pFile[1],"\n          Range            Cycle       Ave     Max     Ave     Min      Max");
	fprintf(pFile[1],"\n         (units)           Counts      Amp     Amp     Mean    Valley   Peak \n");

	for(i=13;i>=1;i--)
    {
    	fprintf(pFile[1],"  %8.2lf to %8.2lf\t%8.1lf\t%6.4g\t%6.4g\t%6.4g\t %6.4g\t %6.4g\n",L[i],L[i+1],C[i],AverageAmp[i],MaxAmp[i],AverageMean[i],MinValley[i],MaxPeak[i]);

     	fprintf(pFile[2]," %ld \t %g \n",i,C[i]);
    }

	fclose(pFile[0]);
	fclose(pFile[1]);
	fclose(pFile[2]);

	printf("\n\n  Total Cycles = %g  hold=%ld  NP=%ld ymax=%g\n",sum,hold,NP,ymax);
	fprintf(pFile[1],"\n\n  Total Cycles = %g  hold=%ld  NP=%ld ymax=%g\n",sum,hold,NP,ymax);

	printf("\n\n The output files are: \n");

    printf(" %s \n",filename[1]);
	printf(" %s \n",filename[2]);

	// printf("\n The dimensions are: Range Units & Cycle Counts \n");

	// printf("\n Please call the output file into your own plotting program. \n");

	printf("\n Press any key to exit.\n");

	getch();

	exit(1);
}