%
% This program solves hover performance of a two-bladed rotor
% using Landgrebe's Prescribed Wake Model
%
b=2;
AR=13.7;
theta1=0.0;
theta0=5.0;
alpha0=0.;
ainf=5.7;
d0=0.011;
d1=0.0;
d2=0.4;
pi=3.14158;
theta0=theta0 *pi/180.;
theta1=theta1*pi/180.;
alpha0=alpha0*pi/180.;
c=1./AR;
rc=0.2*c;
sigma=b*c/pi;
M=90;
n1=36;
n2=20;
deltasi=0.5*pi/M;
si=0;
rn0=1.0;
r=zeros(M);
cl=zeros(M);
v=zeros(M);
dr=zeros(M);
for i=1:M
   si=si+deltasi;
   rn=cos(si);
   dr(i)=rn0-rn;
   r(i)=rn+0.5*dr(i);
   rn0=rn;
end
gamam=0.;
ct=0.;
for i=1:M
   vinf=r(i);
   theta=theta0+r(i)*theta1;
   vi=sigma*ainf*(sqrt(1.+32.*theta*r(i)/sigma/ainf)-1.)/16.;
   alpha = theta - alpha0 - vi/vinf;
   cl(i)=ainf*alpha;
   gama(i)=0.5*vinf*c*cl(i);
   if(gama(i)> gamam)
      gamam=gama(i);
   end
   ct=ct+0.5*sigma*vinf*vinf*cl(i)*dr(i);
end
ct
gamam
r0=0.2;
%
% iterate on CT a few times
%
for itglobal=1:5
   BTIP=1.-(2.*ct)^0.5/b;
   a=0.78;
   lambda=0.145+27.*ct;
   k1=-.25*(ct/sigma);
   k2=-sqrt(ct);
   n=n1*n2;
   dpsi=2.*pi/n1;
   ddpsi=2.*pi/b;
   psi=0;
   psi0=0;
   rv0=1.0;
   zv0=0.0;
   
   for i=1:M
      v(i)=0.0;
   end
   for j=1:n
      psi=j*dpsi;
      rv=a+(1-a)*exp(-lambda*psi);
      if(psi<2.*pi/b)
         zv=k1*psi;
      else
         zv=k1 * 2.*pi/b+ k2 *(psi-2.*pi/b);
      end
      for k=1:b
         psi1 = j*dpsi + (k-1)* 2*pi/b;
         psi0 = psi1 - dpsi;
         xv0=rv0*cos(psi0);
         yv0 = - rv0*sin(psi0);
         xv = rv * cos(psi1);
         yv =- rv * sin(psi1);
         for i = 1:M
            
           
            x0 = r(i);
            x1=xv0;
            y1=yv0;
            z1=zv0;
            x2=xv;
            y2=yv;
            z2=zv;
            r1cr2z=y1*(x0-x2)-y2*(x0-x1);
            r12=(x0-x1)^2+y1*y1+z1*z1;
            r22=(x0-x2)^2+y2*y2+z2*z2;
            r1=(r12)^0.5;
            r2=(r22)^0.5;
            r1r2=r1*r2;
            r1pr2=(x0-x1)*(x0-x2)+y1*y2+z1*z2;
            coef1=(r1+r2)*(1.-r1pr2/r1r2);
            coef2=r12*r22-r1pr2^2+rc*rc*(r12+r22-2.*r1pr2);
      
            coef=coef1/coef2;
            w=r1cr2z*coef;
            v(i)=v(i)+w;
         end
         rv0=rv;
         zv0=zv;
         
      end
   end
   vi0=0.0;
   plot(r,v);
   % iterate on gama
   for itn=1:5
      for i=1:M
         vi = -0.25 * gamam * v(i)/pi;
         vinf=r(i);
         theta=theta0+r(i)*theta1;
         alpha=theta-alpha0-vi/vinf;
         cl(i)=ainf*alpha;
         gama(i)=0.5*vinf*c*cl(i);

         if(r(i)>a)
            if(vi > vi0)
               gaman=gama(i);
            end
         end       
         

         vi0=vi;
      end
      gaman
      err_in_gama=gaman-gamam
      gamam=0.5*(gamam+gaman);
   end


   ct0=ct;
   ct=0.0;
   cq=0.0;
   for i=1:M
      v(i)=-0.25 * gaman * v(i)/pi;
         vinf=r(i);
         theta=theta0+r(i)*theta1;
         alpha=theta-alpha0-vi/vinf;
         cl(i)=ainf*alpha;
         cd=d0+d1*alpha+d2*alpha^2;
         cq=cq+.5*sigma*vinf*(v(i)*cl(i)+vinf*cd)*r(i)*dr(i);
         if(r(i) < BTIP)
            if(r(i)>r0)
               ct=ct+0.5*sigma*vinf*vinf*cl(i)*dr(i);
               end
         end
      end
      ct
      err_ct=ct-ct0;
      ct=0.5*(ct+ct0);
end
%   plot(r,v);
 %  plot(r,gama);
   ct
   cq