clc
clear all
rho=1;
G=6.672;
zt=2;
zb=12;
h=zb-zt;
x0=0;
y0=0;
R=2.5;
dz=.01;
dt=1;
ztr=(zt:dz:zb);
nztr=length(ztr);
te=(0:dt:360);
nte=length(te);
matx=(-30:30);
nx=length(matx);
maty=(-30:30);
ny=length(maty);
xtr=(1:nte);
ytr=(1:nte);
Xtr=(1:nte);
Ytr=(1:nte);
v=(1:nztr);
dg=[1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx 1:nx];
for i=1:nx
for j=1:ny
for k=1:nztr
temp=ztr(k)-zt;
r=(R*temp)/h;
u=0;
if r==0
u=0;
else
for l=1:nte
x=matx(i);
y=maty(j);
t=(-1);
temp=te(l)*(pi/180);
xtr(l)=x0+(r*cos(t*temp));
Xtr(l)=xtr(l)-x;
ytr(l)=y0+(r*sin(t*temp));
Ytr(l)=ytr(l)-y;
end
for l=1:nte-1
xx1=Xtr(l);
xx2=Xtr(l+1);
yy1=Ytr(l);
yy2=Ytr(l+1);
dx=xx1-xx2;
dy=yy1-yy2;
r1=sqrt(xx1^2+yy1^2);
r2=sqrt(xx2^2+yy2^2);
r12=sqrt(dx^2+dy^2);
if r1==0
r1=1e-16;
end
if r2==0
r2=1e-16;
end
if r12==0
r12=1e-16;
end
c=((xx1*xx2)+(yy1*yy2))/(r1*r2);
p=((dy*xx1)+(dx*yy1))/r12;
q=((dx*xx1)+(dy*yy1))/(r1*r12);
f=((dx*xx2)+(dy*yy2))/(r2*r12);
m=((yy1*xx2)-(xx1*yy2))/(r1*r2);
if p>0
s=1;
else
s=-1;
end
if m>0
w=1;
else
w=-1;
end
sd=sqrt(p^2+ztr(k)^2);
if sd==0
sd=1e-16;
end
temp=(ztr(k)*q*s)/sd;
temp1=(ztr(k)*f*s)/sd;
atemp=asin(temp);
atemp1=asin(temp1);
ac=acos(c);
u=u+(w*ac)-atemp-atemp1;
end
temp2=0;
v(k)=G*rho*u;
temp2=temp2+v(k);
end
end
dg(i,j)=((temp2*2)-v(1)-v(nztr))*(dz/2);
end
end
plot(matx,dg(:,31))