急救～牛逼人士进来看下～

nima

#include "stdio.h"
#include "math.h"
#include "iostream.h"
#define g  9.8
#define PI 3.141592
void main()
{ int i,j,N=5,NODES,INDEX;
   float WTT,L=5000.0,D=30.00,F=0.0200,A=4060,AK,VZERO=5.00,VL,VR,HZERO=715.5,HL,HR,Z1=100.0,Z2=50.0,T,TMAX=4.0,TCLOSE=1.0,DELL,DELT,DTDS,C,DELEL,DELHF,SINE,HEADMX,HEADMN;
   float X[20],V[20],H[20],HLOW[20],HHIGH[20],HEAD[20],PIPEZ[20],VNEW[20],HNEW[20];
     WTT=L/A;
     T=0.0;
      C=g/A;
      DELEL=(Z1-Z2)/N;
       DELL=L/N;
       DELT=DELL/(A+VZERO);
   DTDS=DELT/DELL;
      SINE=DELEL/DELL;
      INDEX=TMAX/DELT+1;
       NODES=N+1;
     printf("**water hammer in a simple pipe**\n");
      printf("ALL INPUT DATAS:\n");
      printf(" L/A=%f\n  DELT=%f\n",WTT,DELT);
      printf(" N=5 L=5000.0FT D=30.00IN F=0.0200 A=4060FPS\n",N,L,D,F,A);
      printf("VZERO=5.00FPS HZERO=715.5FT Z1=100.0FT Z2=50.0FT TMAX=4.0 TCLOSE=1.0\n  ",VZERO,HZERO,Z1,Z2,TMAX,TCLOSE);
     AK=(1/C)*F*DELT/(2*D);
     DELHF=12*DELL*VZERO*VZERO/(64.4*D);
    for(i=1;i<=NODES;i++)
      { V=VZERO;
   H=HZERO-(i-1)*DELHF;
    HLOW=H;
     X=(i-1)*DELL/L;
      PIPEZ=Z1+(i-1)*DELEL;
       HEAD=H-PIPEZ;
    }
       printf("**pressure heads,h-values and velocities as functions of time.**\n");
        printf("TIME=%3.3f     X=%3.3f   HEAD-FT=%3.3f   H-FT=%3.3f   V-FPS=%3.3f\n",T,X,HEAD,H,V,i=1,NODES);
  for(j=1;j<=INDEX;j++)
     {  T=T+DELT;
        for(i=2;i<=N;i++)
         { VL=V+DTDS*A*(V[i-1]-V);
        VR=V+DTDS*A*(V[i+1]-V);
        HL=H+DTDS*A*(H[i-1]-H)*(A+VL);
        HR=H+DTDS*A*(H[i+1]-H)*(A-VR);
        VNEW=0.5*(VL+VR+C*(HL-HR)+C*DELT*(VL-VR)*SINE-AK*(VL*abs(VL)+VR*abs(VR)));
        HNEW=0.5*(HL+HR+(VL-VR)/C+DELT*SINE*(VL+VR)-(AK/C)*(VL*abs(VL))-VR*abs(VR));
     }
        HNEW[1]=H[1];
                  VR=V[1]+DTDS*A*(V[2]-V[1]);
          HR=H[1]+DTDS*(H[2]-H[1])*(A-VR);
           VNEW[1]=VR+C*(HNEW[1]-HR)-AK*VR*abs(VR)-C*DELT*VR*SINE;
     if(T>TCLOSE)
     VNEW[NODES]=0.0;
       else
    {VNEW[NODES]=VZERO*(1.0-T/TCLOSE);
      VL=V[NODES]+DTDS*A*(V[N]-V[NODES]);
       HL=H[NODES]+DTDS*(H[N]-H[NODES]*(A+VL));
       HNEW[NODES]=H[N]+(V[N]-VNEW[NODES]-AK*V[N]*abs(V[N]))/C+DELT*VL*SINE;
    }
        for(i=1;i<=NODES;i++)
       { if(HNEW<HLOW)
     HLOW=HNEW;
      if(HNEW>HHIGH)
        HHIGH=HNEW;
      HEAD=HNEW-PIPEZ;
     }
     if(j%i==0)
       printf("TIME=%3.3f\n     X=%3.3f   HEAD-FT=%3.3f   H-FT=%3.3f   V-FPS=%3.3f\n",T,X,HEAD,HNEW,VNEW,i=1,NODES);
      if(T>TMAX)
     printf("table of extreme values\n");
      for(i=1;i<=NODES;i++)
       { V=VNEW;
          H=HNEW;
        }
     
     }
        printf("X            HEADMX           HEADMN           MAX H           MIN H\n");        
        for(i=1;i<=NODES;i++)
       { if(HNEW<HLOW)
     HLOW=HNEW;
     else
     if(HNEW>HHIGH)
     { HHIGH=HNEW;}
        HEAD=HNEW-PIPEZ;
        }
      for(i=1;i<=NODES;i++)
        { HEADMX=HHIGH-PIPEZ;
    HEADMN=HLOW-PIPEZ;
     printf("X=%3.3f   HEADMX=%3.3f   HEADMN=%3.3f  MAX H=%3.3f  MIN H=%3.3f\n",X,HEADMX,HEADMN,HHIGH,HLOW);
         }
        
  }
运行结果总是出错，VNEW出错出的比较离谱