ФЕДЕРАЛЬНАЯ СЛУЖБА ПО ГИДРОМЕТЕОРОЛОГИИ И МОНИТОРИНГУ ОКРУЖАЮЩЕЙ СРЕДЫ (Росгидромет)

//...................................................................................................................................................

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

}