Уважаемые эксперты. Прошу помочь в визуализации. Хочу получить интерактивное динамическое изображение по времени, в коде цикл по времени lv. чуток знаком с gnuplot и файл P.dat для gnuplot. И в изображении я вижу не совсем то что я ожидал там 7 скважин а по задаче 9 должно.

Приложение:
#include<stdio.h>
#include<time.h>
#include<math.h>
#include <cstdlib>
#include <iostream>
#define ronefti 1e3
#define N1 49
#define N2 49
#define N3 20

#define R(i,j,i3) r[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define V(i,j,i3) v[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define W(i,j,i3) w[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define Y(i,j,i3) y[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define YY(i,j,i3) yy[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define F(i,j,i3) f[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define MUX(i,j,i3) mux[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define MUY(i,j,i3) muy[i+(N1+1)*j+(N1+1)*(N2+1)*i3]
#define MUZ(i,j,i3) muz[i+(N1+1)*j+(N1+1)*(N2+1)*i3]

#define NITER 20
#define L1 40000.0
#define L2 40000.0
#define hpl 200.0
#define pi 3.1415926535
#define p0 1.0e7
#define rc 0.1
#define eps 1e-6
#define s0 0.3
#define s1 0.2
#define s2 0.85
#define mu1 1.0e-3
#define mu2 1.0e-2
#define k 1.0e-12
#define m 0.2
#define k1 5
#define k2 4
typedef double type;
static double muxp[N1+1][N2+1][N3+1],muyp[N1+1][N2+1][N3+1];
static double s[N1+1][N2+1][N3+1],ss[N1+1][N2+1][N3+1];
void CreateBuffer(type**,int);

static double qq[k1+k2+1];
static int kx[k1+k2+1],ky[k1+k2+1],lv;
double l;
type *r,*v,*w,*y,*yy,*f,*mux,*muy,*muz;
FILE *tecT;
char *stroka, symbol;
int kn1,kn2,kn3,kn4,propus;
double mu0;

double z1(double x,double y){ return x*(L1-x)*y*(L2-y)*16/L1/L2/L1/L2*hpl/2;}
double z2(double x,double y){ return x*(L1-x)*y*(L2-y)*16/L1/L2/L1/L2*hpl/2+hpl/2;}
double f1(double s)
{ if (s>s1) return exp(3.5*log((s-s1)/0.8)); else return 0.0; }

double f2(double s)
{ if (s<s2) return (1+2.4*s)*exp(2.8*log((s2-s)/s2)); else return 0.0; }

double fi(double s) { double a; a=f1(s); return a/(a+mu0*f2(s)); }

double Q(double s) { return f1(s)/mu1+f2(s)/mu2; }

double epsFikt(double x,double y,double z){if (z>=z1(x,y)&&z<=z2(x,y))return 1; else return 1e-7;}

int main(int argc, char**argv)
{
double e, h, hh,h3,ohh3, a=0, b=0, c, d, tau, tau1, alpha, alpha1, q0;
int i,j,j1,j2,j3,j4,j5,j6,l,i3,n,g;
double t, br, koef, pfix, pfiy, pfiz;
time_t konechvrem,nachvrem; clock_t start,end;
start=clock(); time(&nachvrem);
kn1=0; kn2=N1; kn3=0; kn4=N2; propus=0;

for (l=1; l<=k1; l++) qq[l]=0.8e-5;
for (l=k1+1; l<=k1+k2; l++) qq[l]=-1.0e-5;

kx[1]=N1/2; ky[1]=N1/2;
kx[2]=N1/3; ky[2]=N1/3;
kx[3]=N1/3; ky[3]=2*N1/3;
kx[4]=2*N1/3; ky[4]=N1/3;
kx[5]=2*N1/3; ky[5]=2*N1/3;
kx[6]=N1/2; ky[6]=N1/4;
kx[7]=N1/4; ky[7]=N1/2;
kx[8]=3*N1/4; ky[8]=N1/2;
kx[9]=N1/2; ky[9]=3*N1/4;

h=L1/N1; hh=h*h; h3=hpl/N3;
ohh3=hh/h3/h3;
CreateBuffer(&r,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&v,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&w,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&y,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&yy,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&f,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&mux,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&muy,(N1+1)*(N2+1)*(N3+1));
CreateBuffer(&muz,(N1+1)*(N2+1)*(N3+1));

t=4*m*h3/(1.0e-4);
mu0=mu1/mu2;
// i3=3;
for (i3=0; i3<=N3; i3++) for (i=0; i<=N1; i++) for(j=0; j<=N2; j++)
{ F(i,j,i3)=0.0;R(i,j,i3)=0.0;W(i,j,i3)=0.0;V(i,j,i3)=0.0;
s[i][j][i3]=ss[i][j][i3]=s0; Y(i,j,i3)=p0-h3*i3/9.8*ronefti;YY(i,j,i3)=p0-h3*i3/9.8*ronefti; }
for (lv=1, br=0; lv<30; lv++)
// vot cikl po vremeni
{
t=t*1.1;if (t>1e6) t=1e6;
br+=t; tau=1.0; alpha=1.0; n=0;

koef=k*t/(2*m*hh);
// nachalnye pribligenija
for (i3=0; i3<=N3; i3++)
for (i=0; i<=N1; i++) for(j=0; j<=N2; j++)
{ F(i,j,i3)=0.0;R(i,j,i3)=0.0;W(i,j,i3)=0.0;V(i,j,i3)=0.0;
}

// vychislenie koeff uravnenija filtr i poprav koeff

a=0.5*pi/log(h/rc); b=log(2.5)/log(2.0);
for (i3=0; i3<=N3; i3++) for (i=0; i<=N1; i++) for(j=0; j<=N2; j++)
{muxp[i][j][i3]=k*epsFikt(i*h,j*h,i3*h3);
muyp[i][j][i3]=k*epsFikt(i*h,j*h,i3*h3);}
for (i3=0; i3<=N3; i3++)
for (l=1; l<=k1+k2; l++)
{
i=kx[l]; j=ky[l]; F(i,j,i3)=qq[l]*epsFikt(i*h,j*h,i3*h3);
muxp[i-1][j][i3]=muxp[i-1][j][i3]*a; muxp[i][j][i3]=muxp[i][j][i3]*a;
muxp[i-1][j-1][i3]=muxp[i-1][j-1][i3]*b; muxp[i][j-1][i3]=muxp[i][j-1][i3]*b;
muxp[i-1][j+1][i3]=muxp[i-1][j+1][i3]*b; muxp[i][j+1][i3]=muxp[i][j+1][i3]*b;

muyp[i][j-1][i3]=muyp[i][j-1][i3]*a; muyp[i][j][i3]=muyp[i][j][i3]*a;
muyp[i-1][j-1][i3]=muyp[i-1][j-1][i3]*b; muyp[i-1][j][i3]=muyp[i-1][j][i3]*b;
muyp[i+1][j-1][i3]=muyp[i+1][j-1][i3]*b; muyp[i+1][j][i3]=muyp[i+1][j][i3]*b;
}
// peresylki!!muxp and muyp
// Zadacha B-Lever. Raznostnaja sxema javnyj ugolok protiv potoka
// peresylki ss and y

for (i3=1; i3<N3; i3++)
for (i=1; i<N1; i++) for(j=1; j<N2; j++)
{j1=i-1;j2=i+1;j3=j-1;j4=j+1;j5=i3-1;j6=i3+1;
a=fi(ss[i][j][i3])-fi(ss[i-1][j][i3]); b=fi(ss[i+1][j][i3])-fi(ss[i][j][i3]);
if (Y(j1,j,i3)>Y(j2,j,i3))
pfix=a; else pfix=b;
a=fi(ss[i][j][i3])-fi(ss[i][j-1][i3]); b=fi(ss[i][j+1][i3])-fi(ss[i][j][i3]);

a=fi(ss[i][j][i3])-fi(ss[i][j-1][i3]); b=fi(ss[i][j+1][i3])-fi(ss[i][j][i3]);
if (Y(i,j3,i3)>Y(i,j4,i3))
pfiy=a; else pfiy=b;

a=fi(ss[i][j][i3])-fi(ss[i][j][i3-1]); b=fi(ss[i][j][i3+1])-fi(ss[i][j][i3]);
if (Y(i,j,j5)>Y(i,j,j6))
pfiz=a; else pfiz=b;
s[i][j][i3]=ss[i][j][i3]+koef*Q(ss[i][j][i3])*((Y(j2,j,i3)-Y(j1,j,i3))*pfix
+(Y(i,j4,i3)-Y(i,j3,i3))*pfiy

+ohh3*(Y(i,j,j6)-Y(i,j,j5)+h3*9.8*ronefti)*pfiz
)*epsFikt(i*h,j*h,i3*h3);

};/*
double a1; for(i=0; i<=N1; i++)
{
for(j=0; j<=N2; j++){a1=(s[i][j][1]-0.3)*10000000000;printf("%-7.4f",a1);} printf(" \n");
printf(" \n");};
printf(" \n");*/
for (i3=0; i3<=N3; i3++)
for (n=1; n<=k1; n++)
{ i=kx[n];j=ky[n];
s[i][j][i3]=s[i][j][i3]+t*qq[n]/(m*hh)*epsFikt(i*h,j*h,i3*h3);
if (s[i][j][i3]>s2) s[i][j][i3]=s2*epsFikt(i*h,j*h,i3*h3);
}
//peresylka s zdes
// double a1; for(i=0; i<=N1; i++)
// {
// for(j=0; j<=N2; j++) {a1=(s[i][j][1]-0.3)*100000000000;printf("%-7.4f",a1);} printf(" \n");
//}; printf(" \n");

//vychislenie davlenija p[i][j][i3]
for (i3=0; i3<N3; i3++)
for (i=0; i<N1; i++) for(j=0; j<N2; j++)
{
MUX(i,j,i3)=Q((s[i][j][i3]+s[i+1][j][i3])*0.5)*muxp[i][j][i3];
MUY(i,j,i3)=Q((s[i][j][i3]+s[i][j+1][i3])*0.5)*muyp[i][j][i3];
MUZ(i,j,i3)=Q((s[i][j][i3]+s[i][j][i3+1])*0.5)*k*ohh3*epsFikt(i*h,j*h,i3*h3);
R(i,j,i3)=0.0;W(i,j,i3)=0.0;V(i,j,i3)=0.0;
};
// double a1; printf("%i",lv); printf(" \n");
// for(i=0; i<=N1; i++)
// {
// for(j=0; j<=N2; j++) {a1=(MUX(i,j,4))*100000000;printf("%-7.5f",a1);} printf(" \n");
// printf(" \n");};

//peresylka MUX MUY MUZ

// Iter process sopriagennyx gradientov, predobuslavlivatel Jakobi

n=0; do
{
// vychislenie neviazki i popravki
for (i3=1; i3<N3; i3++) for(i=1; i<N1; i++) for(j=1; j<N2; j++)
{ j1=i-1;j2=i+1;j3=j-1;j4=j+1;j5=i3-1;j6=i3+1;
R(i,j,i3)=MUX(j1,j,i3)*(Y(i,j,i3)-Y(j1,j,i3))+MUX(i,j,i3)*(Y(i,j,i3)-Y(j2,j,i3)) +MUY(i,j3,i3)*(Y(i,j,i3)-Y(i,j3,i3))+MUY(i,j,i3)*(Y(i,j,i3)-Y(i,j4,i3))+ MUZ(i,j,j5)*(Y(i,j,i3)-Y(i,j,j5))+MUZ(i,j,i3)*(Y(i,j,i3)-Y(i,j,j6)) -F(i,j,i3);
W(i,j,i3)=R(i,j,i3)/(MUX(j1,j,i3)+MUX(i,j,i3)+MUY(i,j3,i3)+MUY(i,j,i3)
+MUZ(i,j,j5)+MUZ(i,j,i3)
);
}
for (i=1; i<N1; i++)
for(j=1; j<N2; j++)
{ j1=N3-1;
R(i,j,0)=MUZ(i,j,0)*(Y(i,j,1)-Y(i,j,0));
W(i,j,0)=R(i,j,0)/MUZ(i,j,0);
R(i,j,N3)=MUZ(i,j,j1)*(Y(i,j,j1)-Y(i,j,N3));
W(i,j,N3)=R(i,j,N3)/MUZ(i,j,j1);
}
// Vychislenie v=Aw
for (i3=1; i3<N3; i3++) for(i=1; i<N1; i++) for(j=1; j<N2; j++){ j1=i-1;j2=i+1;j3=j-1;j4=j+1;j5=i3-1;j6=i3+1;
V(i,j,i3)=MUX(j1,j,i3)*(W(i,j,i3)-W(j1,j,i3))+MUX(i,j,i3)*(W(i,j,i3)-W(j2,j,i3))
+MUY(i,j3,i3)*(W(i,j,i3)-W(i,j3,i3))+MUY(i,j,i3)*(W(i,j,i3)-W(i,j4,i3))
+MUZ(i,j,j5)*(W(i,j,i3)-W(i,j,j5))+MUZ(i,j,i3)*(W(i,j,i3)-W(i,j,j6))
;}
for (i=1; i<N1; i++) for(j=1; j<N2; j++)
{j1=N3-1;
V(i,j,0)=MUZ(i,j,0)*(W(i,j,1)-W(i,j,0));
V(i,j,N3)=MUZ(i,j,j1)*(W(i,j,j1)-W(i,j,N3));
}
// Skaliarnoe proizvedenie (w,r) (w,Aw)
b=a; a=0.0;c=0.0;
for (i3=0; i3<=N3; i3++)
for(i=1; i<N1; i++) for(j=1; j<N2; j++) {a=a+W(i,j,i3)*R(i,j,i3);
c=c+W(i,j,i3)*V(i,j,i3);
}
// Vychislenie parametrov iterac processa
tau1=tau; tau=a/c; alpha1=alpha;
if(n>0) alpha=tau1*alpha1*b/(tau1*alpha1*b-tau*a);
// Vychislenie ocherednogo pribligenia
for (i3=0; i3<=N3; i3++)
for(i=1; i<N1; i++) for(j=1; j<N2; j++)
{
e=YY(i,j,i3); YY(i,j,i3)=Y(i,j,i3);
Y(i,j,i3)=alpha*Y(i,j,i3)+(1.0-alpha)*e-alpha*tau*W(i,j,i3);
}

// max|y-yy|
d=0.0;
for (i3=1; i3<N3; i3++)
for (i=1; i<N1; i++) for (j=1; j<N2; j++)
if(fabs(Y(i,j,i3)-YY(i,j,i3))>d) d=fabs(Y(i,j,i3)-YY(i,j,i3));
d=d/p0; n=n+1;
}
// while((n<=NITER)&&(d>=eps));
while(n<NITER);
// vyxod iz iterac cikla
for (i3=0; i3<=N3; i3++) for (i=0; i<=N1; i++) for (j=0; j<=N2; j++) ss[i][j][i3]=s[i][j][i3];
// printf("l=%i n=%i s[i0][j0]=%g br=%g\n",lv,n,s[kx[1]][ky[1]][i3],br/(365*86400));
}

// printf("kol tochek X %d kol tochek Y %d s shagom %f \n", N1, N2, h);
time(&konechvrem); end=clock();
// printf("nachalo scheta: %s \n",ctime(&konechvrem));
// printf("konec scheta: %s \n",ctime(&nachvrem));



// Pechat raspred vodonasyshennosti
/* for(i=0; i<=N1; i++)
{
for(j=0; j<=N2; j++) printf("%-6.3f",s[i][j][1]); printf(" \n");
} printf(" \n"); g=N2/10;
*/
// Pechat raspred davlenia
double a1; for(i=0; i<=N1; i++)
{
for(j=0; j<=N2; j++) { a1=(Y(i,j,2))/p0; printf("%-7.4f",a1); }

printf(" \n");
}
printf("%f %f %d\n",a,c,n);
tecT = fopen("P.dat", "wt");
for(i=0; i<=N1; i++)
{
for(j=0; j<=N2; j++) { a1=(Y(i,j,2))/p0; fprintf(tecT,"%-7.4f",a1); }

fprintf(tecT," \n");
} fclose(tecT);

// i3=N3/2;
system("PAUSE"); return 0;
}
void CreateBuffer(type **buffer,int size)
{
if ((*buffer = (type*)malloc(sizeof(type)*size)) == NULL)
{
printf("\nNot enough memory to allocate buffer!\n");
// MPI_Finalize();
exit(0);
}
}