Алгоритмическая и программная реализация прогноза почвенных параметров

DO 160 J=1,K

JJ=JJ+J

150 SQ=SQRT(R(JJ))

DO 160 I=1,M

L=L+1

160 V(L)=SQ*V(L)

RETURN

END

SUBROUTINE TRACE (M,R,CON,K,D)

DIMENSION R(1),D(1)

FM=M

L=0

DO 100 I=1,M

L=L+I

cwrite(6,211) r(l),con

211format(1x,'tr r=',f12.5,' con=',f12.5)

cpause

100 D(I)=R(L)

K=0

DO 110 I=1,M

IF(D(I)-CON) 120,105,105

105 K=K+1

110 D(I)=D(I)/FM

120 DO 130 I=2,K

130 D(I)=D(I)+D(I-1)

cwrite(6,212) k

cpause

212format(1x,'tr w k=',i5)

RETURN

END

SUBROUTINE VARMX (M,K,A,NC,TV,H,F,D)

DIMENSION A(1),TV(1),H(1),F(1),D(1)

EPS=0.00116

TVLT=0.0

LL=K-1

NV=1

NC=0

FN=M

FFN=FN*FN

CONS=0.7071066

DO 110 I=1,M

H(I)=0.0DO 110 J=1,K

L=M*(J-1)+I

110 h(i)=h(i)+a(l)*a(l)

DO 120 I=1,M

115 h(i)=sqrt(h(i))

DO 120 J=1,K

L=M*(J-1)+I

120 A(L)=A(L)/H(I)

GO TO 132

130 NV=NV+1

TVLT=TV(NV-1)

132 TV(NV)=0.0

DO 150 J=1,K

AA=0.0

BB=0.0

LB=M*(J-1)

DO 140 I=1,M

L=LB+I

CC=A(L)*A(L)

AA=AA+CC

140 bb=bb+cc*cc

150 TV(NV)=TV(NV)+(FN*BB-AA*AA)/FFN

IF(NV-51) 160, 430, 430160 IF((TV(NV)-TVLT)-(1.E-7)) 170, 170, 190

170 NC=NC+1

IF(NC-3) 190, 190, 430

190 DO 420 J=1,LL

L1=M*(J-1)

II=J+1

DO 420 K1=II,K

L2=M*(K1-1)

AA=0.0

BB=0.0

CC=0.0

DD=0.0

DO 230 I=1,M

L3=L1+I

L4=L2+I

U=(A(L3)+A(L4))*(A(L3)-A(L4))

T=A(L3)*A(L4)

T=T+T

CC=CC+(U+T)*(U-T)

DD=DD+2.0*U*T

AA=AA+U

230 BB=BB+T

T=DD-2.0*AA*BB/FN

B=CC-(AA*AA-BB*BB)/FN

IF(T-B) 280, 240, 320

240 IF((T+B)-EPS) 420,250,250

250 COS4T=CONS

SIN4T=CONS

GO TO 350

280 TAN4T=ABS(T)/ ABS(B)

IF(TAN4T-EPS) 300,290,290

290 COS4T=1.0/ SQRT(1.0+TAN4T*TAN4T)

SIN4T=TAN4T*COS4T

GO TO 350

300 IF(B) 310,420,420

310 SINP=CONS

COSP=CONS

GO TO 400

320 CTN4T= ABS(T/B)

IF(CTN4T-EPS) 340,330,330

330 SIN4T=1.0/ SQRT(1.0+CTN4T*CTN4T)

COS4T=CTN4T*SIN4T

GO TO 350

340 COS4T=0.0

SIN4T=1.0

350 COS2T=SQRT((1.0+COS4T)/2.0)

SIN2T=SIN4T/(2.0*COS2T)

355 COST=SQRT((1.0+COS2T)/2.0)

SINT=SIN2T/(2.0*COST)

IF(B) 370,370,360

360 COSP=COST

SINP=SINT

GO TO 380

370 COSP=CONS*COST+CONS*SINT

375 SINP=ABS(CONS*COST-CONS*SINT)

380 IF(T) 390,390,400

390 SINP=-SINP

400 DO 410 I=1,M

L3=L1+I

L4=L2+I

AA=A(L3)*COSP+A(L4)*SINP

A(L4)=-A(L3)*SINP+A(L4)*COSP

410 A(L3)=AA

420 CONTINUE

GO TO 130

430 DO 440 I=1,M

DO 440 J=1,K

L=M*(J-1)+I

440 A(L)=A(L)*H(I)

NC=NV-1

DO 450 I=1,M

450 H(I)=H(I)*H(I)

DO 470 I=1,M

F(I)=0.0

DO 460 J=1,K

L=M*(J-1)+I

460 F(I)=F(I)+A(L)*A(L)

470 D(I)=H(I)-F(I)

RETURN

END

subroutine grpk(m,k1,v,mnp,ip1,mnpi,kis)

DIMENSION V(1)

dimension mnf(151)

dimension mnp(1),mnpi(1),mrez(151)

dimension zna(151),nf(151),mgp(151)

common anaz(151,60),knaz

cwrite(6,7001) m,kis

7001format(1x,'m=',i3,' kis=',i3)

IKR=0

IF(IP1.NE.0) GO TO 91

DO 92 I=1,KIS

92 MREZ(I)=MNPI(I)

GO TO 93

91 CONTINUE

cwrite(6,4501) (mnp(i),i=1,ip1)

4501format(1x,'mnp',15i4)

DO 1 IS=1,KIS

DO 2 IP=1,IP1

IF(IS.EQ.MNP(IP)) GO TO 1

2 CONTINUE

IKR=IKR+1

MREZ(IKR)=MNPI(IS)

1 CONTINUE

cwrite(6,4502) (mnpi(i),i=1,kis)

4502format(1x,'mnpi',15i4)

cwrite(6,4503) (mrez(i),i=1,ikr)

4503format(1x,'mrez',15i4)

93 CONTINUE

KT=M*K1

DO 5 I=1,M

XMAX=ABS(V(I))

INS=I

DO 6 J=I,KT,M

IF(XMAX-ABS(V(J))) 7,6,6

7 CONTINUE

XMAX=ABS(V(J))

INS=J

6 CONTINUE

JST=INS/M+1

S1=INS*1./M

I1=S1

IF(S1.EQ.I1) JST=I1

ZNA(I)=V(INS)

NF(I)=JST

5 CONTINUE

DO 10 I=1,M

NFT=NF(I)

MNF(I)=NFT

IF(I.EQ.1) GO TO 201

II=I-1

DO 206 IR=1,II

IF(NFT.EQ.MNF(IR)) GO TO 10

206 CONTINUE

201 CONTINUE

IB=0

DO 11 J=I,M

IF(NFT.NE.NF(J)) GO TO 11

IB=IB+1

MGP(IB)=J

11 CONTINUE

cwrite(6,4504) (mgp(iu),iu=1,ib)

4504format(1x,'mgp',15i4)

WRITE(6,12)NFT

12 FORMAT(/9X,'Таблица-Объединение по фактору',I5,'.'/)

write(6,7013)

7013format(5x,58('-'))

WRITE(6,7012)

write(6,7013)

7012format(5x, '|НОМЕР| НАЗВАНИЕ ПАРАМЕТРА |',

*' НАГРУЗКА |')

ir=mgp(1)

ismax=mrez(ir)

xmax=abs(zna(ir))

DO 13 IT=1,IB

IR=MGP(IT)

IS=MREZ(IR)

if(xmax.ge.abs(zna(ir))) goto 12102

ismax=is

xmax=abs(zna(ir))

12102continue

cif(knaz.gt.70) goto 801

do 214 j=35,60

cwrite(*,287) j,anaz(i,j)

287format(1x,'j=',i4,60a1)

if(anaz(is,j).eq.' ') goto 214

cwrite(*,287) j,(anaz(is,js),js=1,60)

cpause

goto 215

214continue

write(6,217) is,(anaz(is,j),j=1,34),zna(ir)

217format(5x,'|',i4,' |',34a1,'|',f11.4,4x,'|')

goto 218

215 WRITE(6,14) is,(anaz(is,j),j=1,60),ZNA(IR)

218write(6,7013)

14 FORMAT(5X,'|',i4,' |',34a1,'| |'/5x,

* '| |',26a1,8x,'|',F11.4,4x,'|')

cgoto 13

c801write(6,141) is,zna(ir)

141 format(10x,'показатель-',i4,' значение нагрузки= ',f10.5)

13 CONTINUE

write(6,12101) nft,ismax,(anaz(ismax,j),j=1,60)

12101format(1x,'в факторе-',i5,' базовый параметр-',i5/1x,60a1/)

10 CONTINUE

RETURN

END

c-умножение матрицы a на матрицу b

subroutine gtprd(a,b,r,n,m,l)

dimension a(1),b(1),r(1)

ir=0

ik=-n

do 10 k=1,l

ij=0

ik=ik+n

do 10 j=1,m

ib=ik

ir=ir+1

r(ir)=0.

do 10 i=1,n

ij=ij+1

ib=ib+1

10 r(ir)=r(ir)+a(ij)*b(ib)

return

end

c-метод Гаусса решения систем

subroutine gay(n,a,b,x)

dimension a(1),b(1),x(1)

cdouble precision a,b,x,t,tz,p,sum

i=1

81j=i

ct=a(i,i)

t=a((i-1)*n+i)

if(t.ne.0) goto 2

l=i+1

4if(a(l+(j-1)*n).ne.0) goto 3

l=l+1

if(l.le.n) goto 4

write(6,*) 'матрица вырожденная'

stop

3continue

do 5 js=1,n

tz=a(l+(js-1)*n)

a(l+(js-1)*n)=a(i+(js-1)*n)

a(i+(js-1)*n)=tz

5continue

t=a(i+(i-1)*n)

2continue

a(i+(j-1)*n)=a(i+(j-1)*n)/t

j=j+1

if(j.le.n) goto 2

b(i)=b(i)/t

k=i+1

71p=a(k+(i-1)*n)

j=i

7a(k+(j-1)*n)=a(k+(j-1)*n)-p*a(i+(j-1)*n)

j=j+1

if(j.le.n) goto 7

b(k)=b(k)-p*b(i)

k=k+1

if(k.le.n) goto 71

i=i+1

if(i.lt.n) goto 81

cwrite(6,*) ' матрица после прямого хода'

do 56 i=1,n

cwrite(6,57) (a(i+(j-1)*n),j=1,n)

57format(1x,5f12.5)

56continue

if(a(n+(n-1)*n).ne.0) goto 82

write(6,*) 'матрица вырожденная'

stop

82continue

nt=n

x(n)=b(n)/a(n+(n-1)*nt)

cwrite(6,86) x(n)

85n=n-1

if(n.lt.1) goto 83

sum=0

j=n+1

84sum=sum+a(n+(j-1)*nt)*x(j)

j=j+1

if(j.le.nt) goto 84

x(n)=b(n)-sum

cwrite(6,86) x(n)

goto 85

83continue

cwrite(5,86)(x(i),i=1,nt)

86format(1x,5f15.5)

return

end

Набор normmatr.cpp(построение нормализованной матрицы)

#include<windows.h>

#include<string.h>

#include<stdio.h>

#include<alloc.h>

#include<dos.h>

#include"menu.h"

#include<math.h>

long int soch()

{

extern FILE *fp3,*fp5,*fp8;

extern void randu(long int ix,long int *piy,float *pfr);

extern void norm(float sr,float si,float fr,float *pzn);

float srm[201],sigm[201],fr,zn,sr,si;

char zgl[81];

long int i,m,kg,ix,iy,ig,j,isc;

//open(3,file='isxnorm',status='old')

//open(5,file='wixnorm',status='Unknown')

fp3=fopen("isxncpp","r");

fp5=fopen("wixncpp","w");

fp8=fopen("slych","w");

for(i=1;i<=80;i++)

{

fscanf(fp3,"%c",&zgl[i]);

fprintf(fp5,"%c",zgl[i]);

}

fscanf(fp3,"%ld%ld",&kg,&m);

fprintf(fp5,"\n%ld %ld\n",kg,m);

for(i=1;i<=m;i++)

{

fscanf(fp3,"%f%f",&srm[i],&sigm[i]);

//fprintf(fp5,"%f %f\n",srm[i],sigm[i]);

}

ix=55555l;

isc=0;

for(ig=1;ig<=kg;ig++)

{

for(j=1;j<=m;j++)

{

randu(ix,&iy,&fr);

fprintf(fp8,"%ld ",iy);

isc=isc+1;

if(isc<5) goto r1;

isc=0;

fprintf(fp8,"\n");

r1:;

ix=iy;

sr=srm[j];

si=sigm[j];

norm(sr,si,fr,&zn);

//rab[j]=zn;

fprintf(fp5,"%f ",zn);

}

fprintf(fp5,"\n");

}

fclose(fp3);

fclose(fp5);

fclose(fp8);

return(1l);

}

void norm(float sr,float si,float alf,float *px)

{

float xg,xn,xn1,s,dx,zk,xt,ar,x;

//c-sr,srednee

//c-si,sqrt(d)

//c-alf,znachenie funkcii raspredelenij

//c-alf rawnomerno raspredeleno

//c-xg,max pokazatel exp

//c-xn,inf integrirowanij

//c-xg=9

//c-s,tekyshee znachenie integrala

xg=18;

xn=(sr-pow(2.*si,0.5)*pow(1.*si*xg,0.5));

xn1=sr+pow(2.*si,0.5)*pow(1.*si*xg,0.5);

s=0.;

//c-xt,argument

//c-dx,shag integrirowanij

//c-ЗДЕСЬ ВСЕ ВЕРНО,УВАЖАЕМЫЕ!

dx=(xn1-xn)/10000;

zk=1./(pow(8.*atan(1.),0.5)*si);

//cwrite(*,651) sr,si,xn1,xn,dx,zk

//651format(1x,'norm','sr=',e15.7,' si=',e15.7,' xn1=',e15.7/1x,

// *' xn=',e15.7,' dx=',e15.7,' zk=',e15.7)

xt=xn;

m1: ar=(xt-sr)/si*(xt-sr)/(2*si);

//c-ar,argument exp

s=s+exp(-ar)*dx*zk;

//c-zk,koeff

if(s>=alf) goto m2;

xt=xt+dx;

if(xt<xn1) goto m1;

m2:;

x=xt;

*px=x;

return;

}

void randu(long int ix,long int *piy,float *pyfl)

{

long int i1,i2,iz,iy,j;

float yfl;

i1=ix/32700l;

i2=ix-i1*32700l;

if(i1==0) goto m25;

iz=32700l*65539l;

iy=0;

for(j=1;j<=i1;j++)

iy=iy+iz;

iy=iy+i2*65539l;

goto m26;

m25: iy=ix*65539l;

m26:;

if(iy<0) goto m5;

if(iy>=0) goto m6;

m5: iy=iy+2147483647l+1l;

m6: yfl=iy;

yfl=yfl*0.4656613e-9;

*piy=iy;

*pyfl=yfl;

return;

}

Набор menu.h

#define IDM_RUN 108

Набор menu.rc

;EXAMPLE FILE RECOURSE FOR MENU

#include<windows.h>

#include"menu.h"

MYMENU MENU

{

MENUITEM "run",IDM_RUN

}

Размещено на Allbest.ru