//...................................................................................................................................................
#include <stdio.h>
#include <math.h>
#define SI 184
#define N 47 //длина исходного ряда
#define SIZ 185
FILE *REZ;
float power(float a, int b);
float a[SI];
float H;
float DEL_DX=1;
float C,S,xx;
float A[SI][SIZ], X[SIZ], YY[N],YYY[N];
//Входные данные. Временной ряд метеопараметров (сумма среднегодовых осадков)
float Y[N]={377.4. 520.4. 445.6. 571.2. 548.1. 599.3. 639. 432. 571. 521. 388. 592. 586. 480. 502. 706.2. 489.9. 489.7. 596.5. 690. 671. 454. 673. 539.9. 427.5. 476. 556.2. 600. 574.3. 659.1. 541. 359.8. 654.4. 642.7. 833.3. 570.7. 560.7. 686. 595.7. 395. 709. 492. 659. 486. 687. 590. 515};
float XX[N]={0.0. 1.0. 2.0. 3.0. 4.0. 5.0. 6.0. 7.0. 8.0. 9.0. 10.0.
11.0.12.0.13.0.14.0.15.0.16.0.17.0.18.0.19.0.20.0.21.0.22.0.23.0.24.0.25.0.
26.0.27.0.28.0.29.0.30.0.31.0.32.0.33.0.34.0.35.0.36.0.37.0.38.0.
39.0.40.0.41.0.42.0.43.0.44.0.45.0.46.0};
int main ()
{
int i, j, k,l, k1,k2,l1,l2,m1, m2, i1, i2, j1, j2, m, n,jj;
H=0.99999/DEL_DX;
for(i=0; i<SI; i++)
{
for(j=0; j<(SI+1); j++)
{
A[i][j]=0.0;
}
}
i=2;
j=2;
k=0;
AZ1: j1=j+1;
A[i][j]=2.0;
A[i][j1]=6*XX[k];
i=i+4;
j=j+4;
k=k+1;
if (i<SI) goto AZ1;
i=6;
j=2;
k=1;
AZ2: j1=j+1;
A[i][j]=-2.0;
A[i][j1]=-6*XX[k];
i=i+4;
j=j+4;
k=k+1;
if (i<SI) goto AZ2;
A[SI-1][SI-2]=2.0;
A[SI-1][SI-1]=6.0*XX[SI/4];
k=0;
l=0;
m=0;
AZ3: i=l;
i1=l+1;
k1=k;
k2=k+4;
m1=m+1;
for(j=k1; j<k2; j++)
{
A[i][j]=power(XX[m].(j-k));
A[i1][j]=power(XX[m1].(j-k));
}
k=k+4;
l=l+4;
m=m+1;
if(k<SI) goto AZ3;
k=0;
l=3;
m=1;
AZ4: i=l;
k1=k+1;
k2=k+4;
for(j=k1; j<k2; j++)
{
j1=j+4;
A[i][j]=(j-k)*power(XX[m].(j-k-1));
A[i][j1]=-(j1-k2)*power(XX[m].(j1-k2-1));
}
k=k+4;
l=l+4;
m=m+1;
if(l<(SI-3)) goto AZ4;
for(i=0.m=0; i<SI; i+=4)
{
i1=i+1;
m1=m+1;
A[i][SI]=Y[m];
A[i1][SI]=Y[m1];
m++;
}
for (i=0; i<(SI-1); i++)
{
C=A[i][i];
for (j=1; j<SIZ; j++)
{
if((j-i) <= 0) continue;
if((j-i) > 0)
*(*(A+i)+j) =(*(*(A+i)+j))/(*(*(A+i)+i));
for (k=1; k<SI; k++)
{
if((k-i)<=0) continue;
if((k-i) > 0)
{
*(*(A+k)+j) = (*(*(A+k)+j)) - (*(*(A+i)+j)) *( *(*(A+k)+i));
}
}
}
}
*(X+(SI-1)) = *(*(A+(SI-1))+SI)/(*(*(A+(SI-1))+(SI-1)));
for(i=0; i<(SI-1); i++)
{
for(j=0. C=0.; j<(SI-1); j++)
{
if((i-j) < 0) continue;
if((i-j) >= 0)
C += *(*(A+(SI-2-i))+(SI-1-j))*(*(X+SI-1-j));
}
*(X+(SI-2-i)) = (*(*(A+(SI-2-i))+(SI))-C);
}
REZ=fopen("rez.dat"."w");
xx=0.0;
for(i=0; i<(N-1); i++)
{
double c1.c2.c3.c4;
n=floor(xx);
c1=X[4*n];
c2=X[4*n+1];
c3=X[4*n+2];
c4=X[4*n+3];
YY[i]=c1+c2*xx+c3*power(xx.2)+c4*power(xx.3);
YYY[i]=c2+2*c3*xx+3*c4*power(xx.2);
//Выходные данные.
//Значения фазового портрета исходного временного ряда
fprintf(rez." %.2f %.fe \n".YY[i].YYY[i]);
xx=xx+H;
}
fclose (REZ);
float power(float a. int b)
{
int i;
float T;
T=1.;
for(i=1; i<=b; i++)
{
T=T*a;
}
return (T);
}
