clear all;
clc;
for t=1:100
 Traffic(t) =curve(t);
end
MaxTraffic = max(Traffic);
w = 0.2;
Wmax = 2*pi*w/3000;
x=[0:10:300];
y=[0:10:300];
Nx=length(x);
Ny=length(y);
Sigma = 0.53;
t = 0；
M = 10*curve(t)/MaxTraffic;
sum = 0;
for i=1:Nx
 forj=1:Ny
  Mu = log(M)-0.5*Sigma^2;
  Rho(i,j) = RhoFromCoordination(x(i),y(j),Wmax,Sigma,Mu);
  Lognrnd(i,j) = round(exp(Sigma*Rho(i,j)+Mu));
  sum = Lognrnd(i,j)+sum;
 end
end
sum
[xi,yi]=meshgrid(0:2:300,0:2:300);
z1=interp2(x,y,Lognrnd,xi,yi,'spline');%三次样条插值
surf(xi,yi,z1)

x=0:100:1200;
y=0:100:1000;
[x y]=meshgrid(x,y);
z=mydata1;
x1=0:10:1200;
y1=0:10:1000;
[x1 y1]=meshgrid(x1,y1);
z1=interp2(x,y,z,x1,y1)
surf(x1,y1,z1)
shading flat
square=0;
[r c]=size(z1)
for n=1:c-1
 for m=1:r-1
  a=x1(m,n+1)-x1(m,n);
  b=y1(m+1,n)-y1(m,n);
  temp=z1(m+1,n+1)-z1(m,n);
  c=sqrt(a*a+b*b+temp*temp);
  temp=z1(m,n+1)-z1(m,n);
  a1=sqrt(temp*temp+a*a);
  temp=z1(m+1,n+1)-z1(m,n+1);
  b1=sqrt(temp*temp+b*b);
  p=(a1+b1+c)/2;
  square=square+sqrt(p*(p-a1)*(p-b1)*(p-c));
  temp=z1(m+1,n)-z1(m,n);
  b2=sqrt(b*b+temp*temp);
  temp=z1(m+1,n+1)-z1(m+1,n);
  a2=sqrt(temp*temp+a*a);
  p=(a2+b2+c)/2;
  square=square+sqrt(p*(p-a2)*(p-b2)*(p-c));
 end
end
square