﻿function nll
format longg
span=2;
ym=-span/2; %%%general lift distribution
yp=span/2;
noseg=101;
dy=(yp-ym)/(noseg-1);
vel=30;
locc=1000000; %%bu neden 1000000
% cdif=1; 
alfa=0.1;  

for i=1:noseg
y(i)=ym+(i-1)*dy;  
c(i)=10*sqrt(1-(2*y(i)/span)^2); %%% c dediğimiz circulation. 10 bizim c0 ımız

end
y';
c';
%%
for i=1:(noseg-1)/2.+1
cbar(i)=0.3+(0.004*(i-1));
%cbar(i)=0.4;
end

for i=1:noseg
cbar(noseg+1-i)=cbar(i);
%cbar(i)=0.1;
end
%%%%%%%%%%%%%
cbar';

iter=0;
while (iter<=200)
iter=iter+1;
% cdif=0;
    
totint=0;
for j=1:noseg
 
 derc=0;
 down=0;
 intpart=0;
 
for i=1:noseg
if (i==1)
derc(i)=(c(i+1)-c(i))/(y(i+1)-y(i));
down(i)=y(j)-y(i);
intpart(i)=derc(i)/down(i);
elseif (i==noseg)
derc(i)=(c(i)-c(i-1))/(y(i)-y(i-1));
down(i)=y(j)-y(i);
intpart(i)=derc(i)/down(i);
    else
derc(i)=(c(i+1)-c(i-1))/(y(i+1)-y(i-1));
down(i)=y(j)-y(i);
intpart(i)=derc(i)/down(i);       
    end        
% locc;
if(abs(down(i))<(dy/10))
locc=i;
end
% locc;

end
% j;
% locc;

if (locc==1)  
intpart(locc)=intpart(locc+1);
elseif(locc==noseg)
    intpart(locc)=intpart(locc-1);
else
        intpart(locc)=(intpart(locc-1)+intpart(locc+1))/2.;
end

% locc=10000000;



totint(j)=0;
for i=2:2:noseg-1;
totint(j)=totint(j)+dy/3*(intpart(i-1)+4*intpart(i)+intpart(i+1));
end

% totint(j);
aind(j)=totint(j)/(4.*acos(-1.)*vel);
aeff(j)=alfa-aind(j);
cl(j)=6.917*(aeff(j));

cnew(j)=0.5*vel*cbar(j)*cl(j);
j;



end
cold=c;
cold';
for i=1:noseg
c(i)=cold(i)+0.05*(cnew(i)-cold(i));
end
c';
aind';
cdif=1;
end

aind';
Lift=0;
indrag=0;
for i=2:2:noseg-1
Lift=Lift+dy*(c(i-1)+4*c(i)+c(i+1))/3.;
% indrag(2)=dy*(c(2-1)*aind(2-1)+4*c(2)*aind(2)+c(2+1)*aind(2+1))/3.;
indrag=indrag+dy*(c(i-1)*aind(i-1)+4*c(i)*aind(i)+c(i+1)*aind(i+1))/3.;
% indrag2(i)=indrag(i)*1.225*vel
end
Lift=Lift*1.225*vel;
indrag=indrag*1.225*vel
S=(cbar(1)+cbar((noseg-1)/2+1))/2.*span;
AR=span^2/S;
CL=Lift*2/(1.225*vel^2*S);
% cdi=indrag*2/(1.225*vel^2*S);
% oswald=CL^2/cdi/acos(-1.)/AR;

end