clear all; 
close all;
clc;
load constsandscaling.mat

for k=1:Nz
  % w = atmoshwm(45,45,1000*Hnorm*z(k),'day',102,'seconds',82800,'apindex',80,'model','total', 'version', '14');
   %w = atmoshwm(45,45,1000*Hnorm*z(k),'day',102,'seconds',82800,'model','disturbance');
    w = atmoshwm([45;45],[45;45],[1000*Hnorm*z(k);1000*Hnorm*z(k)],'day',[102;102],'seconds',[82800-1*3600;82800-1*3600]);
    VnX(k)=3.6*w(1,1)/Velnorm;
    VnY(k)=3.6*w(1,2)/Velnorm;
    VnZ(k)=0;
    Uwind(k)=sqrt(VnX(k)*VnX(k)+VnY(k)*VnY(k));
end
plot(Uwind(kmin:kmax)*Velnorm/3.6, Hnorm*z(kmin:kmax))

 for h=kmin:kmax
      ViZnop(h)=(-VnX(h)*SINI*COSI-VnY(h)*kappaspnormmsise(251,h)*COSI+VnZ(h)*(kappaspnormmsise(251,h)*kappaspnormmsise(251,h)+SINI*SINI))/(1+kappaspnormmsise(251,h)* kappaspnormmsise(251,h));
      Cx(h)=(-SINI*COSI/(1+kappaspnormmsise(251,h)* kappaspnormmsise(251,h)));
      Cy(h)=(-kappaspnormmsise(251,h)*COSI/(1+kappaspnormmsise(251,h)* kappaspnormmsise(251,h)));
    %ViZnop(h)=Cx(h)*VnX(h)+Cy(h)*VnY(h);
 end
    
     for h=kmin+1:kmax-1
  Ceffectivenorm(h)=-ViZnop(h);
  
  Cxprimnorm(h)=(Cx(h+1)-Cx(h))/(z(h+1)-z(h));
  Cyprimnorm(h)=(Cy(h+1)-Cy(h))/(z(h+1)-z(h));
  Cprimeffectivenorm(h)=-VnX(h)*Cxprimnorm(h)-VnY(h)*Cyprimnorm(h);
  divergenceZnorm(h)=(ViZnop(h+1)-ViZnop(h))/(z(h+1)-z(h));
  sheardivergencenorm(h)=Cx(h)*(VnX(h+1)-VnX(h))/(z(h+1)-z(h))+Cy(h)*(VnY(h+1)-VnY(h))/(z(h+1)-z(h));
  shearXnorm(h)=(VnX(h+1)-VnX(h))/(z(h+1)-z(h));
  shearYnorm(h)=(VnY(h+1)-VnY(h))/(z(h+1)-z(h));
    end
    
        for h=kmin+1:kmax-1
  Ceffective(h)=-ViZnop(h)*Velnorm/3.6;
  convergenceZ(h)=- divergenceZnorm(h)/(tnorm*3600);
  shearconvergence(h)=-sheardivergencenorm(h)/(tnorm*3600);
  shearX(h)=shearXnorm(h)/(tnorm*3600);
  shearY(h)=shearYnorm(h)/(tnorm*3600);
  Cprimeffective(h)=Cprimeffectivenorm(h)/(tnorm*3600);
    end
  
    CeffectiveforHczero1=min(abs(Ceffective(550:590)));
    CeffectiveforHczero2=min(abs(Ceffective(590:650)));
    CeffectiveforHczero3=min(abs(Ceffective(650:750)));
  CprimeffectiveforHcprimeffectivemax=max(Cprimeffective);
   convergenceZforHconvergenceZmax=max(convergenceZ);
   shearconvergenceforHshearcovergencemax=max(shearconvergence(h));
    
           for h=550:590
        if abs(Ceffective(h))==CeffectiveforHczero1
            hczero1=h;
                        Hczero1=Hnorm*z(h);
                        directionczero1=(180/pi)*acos((VnX(h)/Uwind(h)));
                        Cprimeffectiveczero1=Cprimeffective(h);
                        convergenceZczero1=convergenceZ(h);
                        shearconvergenceczero1=shearconvergence(h);
                        shearXczero1=shearX(h);
                        shearYczero1=shearY(h);
                        
        end
      end
      
      if VnY(hczero1)<0
          
     directionczero1=360-(180/pi)*acos((VnX(hczero1)/Uwind(hczero1)));
      end
      
      
         for h=590:650
        if abs(Ceffective(h))==CeffectiveforHczero2
                        hczero2=h;
                        Hczero2=Hnorm*z(h);
                        directionczero2=(180/pi)*acos((VnX(h)/Uwind(h)));
                        Cprimeffectiveczero2=Cprimeffective(h);
                        convergenceZczero2=convergenceZ(h);
                        shearconvergenceczero2=shearconvergence(h);
                        shearXczero2=shearX(h);
                        shearYczero2=shearY(h);
                        
        end
         end
         
         if VnY(hczero2)<0
             directionczero2=360-(180/pi)*acos((VnX(hczero2)/Uwind(hczero2)));
         end
      
         for h=650:750
        if abs(Ceffective(h))==CeffectiveforHczero3
                        hczero3=h;
                        Hczero3=Hnorm*z(h);
                        directionczero3=(180/pi)*acos((VnX(h)/Uwind(h)));
                        Cprimeffectiveczero3=Cprimeffective(h);
                        convergenceZczero3=convergenceZ(h);
                        shearconvergenceczero3=shearconvergence(h);
                        shearXczero3=shearX(h);
                        shearYczero3=shearY(h);
                        
        end
      end
        if VnY(hczero3)<0
          directionczero3=360-(180/pi)*acos((VnX(hczero3)/Uwind(hczero3)));
        end
       
       
       
       
      
       for h=kmin+1:kmax-1
        if Cprimeffective(h)==CprimeffectiveforHcprimeffectivemax
            hcprimeffectivemax=h;
            Hcprimeffectivemax=Hnorm*z(h);
            directioncprimeffectivemax=(180/pi)*acos((VnX(h)/Uwind(h)));
                        Cprimeffectivecprimeffectivemax=Cprimeffective(h);
                        convergenceZcprimeffectivemax=convergenceZ(h);
                        shearconvergencecprimeffectivemax=shearconvergence(h);
                        shearXcprimeffectivemax=shearX(h);
                        shearYcprimeffectivemax=shearY(h);
        end
       end  
      
     if VnY(hcprimeffectivemax)<0
           directioncprimeffectivemax=360-(180/pi)*acos((VnX(hcprimeffectivemax)/Uwind(hcprimeffectivemax)));
       end   
       
       
       
       for h=kmin+1:kmax-1
        if shearconvergence(h)== shearconvergenceforHshearcovergencemax
            hshearconvergencemax=h;
            Hshearconvergencemax=Hnorm*z(h);
            directionshearconvergencemax=(180/pi)*acos((VnX(h)/Uwind(h)));
                        Cprimeffectiveshearconvergencemax=Cprimeffective(h);
                        convergenceZshearconvergencemax=convergenceZ(h);
                        shearconvergenceshearconvergencemax=shearconvergence(h);
                        shearXshearconvergencemax=shearX(h);
                        shearYshearconvergencemax=shearY(h);

        end
       end        
       
        if VnY(hshearconvergencemax)<0
           directionshearconvergencemax=360-(180/pi)*acos((VnX(hshearconvergencemax)/Uwind(hshearconvergencemax)));
       end
       
for h=kmin+1:kmax-1
        if convergenceZ(h)== convergenceZforHconvergenceZmax
            hconvergenceZmax=h;
            HconvergenceZmax=Hnorm*z(h);
            directionconvergenceZmax=(180/pi)*acos((VnX(h)/Uwind(h)));
                        Cprimeffectivecprimeffectivemax=Cprimeffective(h);
                        convergenceZconvergenceZmax=convergenceZ(h);
                        shearconvergenceconvergenceZmax=shearconvergence(h);
                        shearXconvergenceZmax=shearX(h);
                        shearYconvergenceZmax=shearY(h);

        end
       end        
      
  if VnY(hconvergenceZmax)<0
     directionconvergenceZmax=360-(180/pi)*acos((VnX(hconvergenceZmax)/Uwind(hconvergenceZmax)));
 end
          disp('Hczero1')
       disp(Hczero1)
       disp('directionczero1')
       disp(directionczero1)
       disp('shearconvergenceczero1')
       disp(shearconvergenceczero1)
       disp('shearXczero1')
       disp(shearXczero1)
       disp('shearYczero1')
       disp(shearYczero1)
       
        disp('Hczero2')
       disp(Hczero2)
       disp('directionczero2')
       disp(directionczero2)
       disp('shearconvergenceczero2')
       disp(shearconvergenceczero2)
       disp('shearXczero2')
       disp(shearXczero2)
       disp('shearYczero2')
       disp(shearYczero2)
       
        disp('Hczero3')
         disp(Hczero3)
       disp('directionczero3')
       disp(directionczero3)
       disp('shearconvergenceczero3')
       disp(shearconvergenceczero3)
       disp('shearXczero3')
       disp(shearXczero3)
       disp('shearYczero3')
       disp(shearYczero3)
       
    
    
       disp('Hcprimeffectivemax')
       disp(Hcprimeffectivemax)
       disp('directioncprimeffectivemax')
       disp(directioncprimeffectivemax)
       disp('shearconvergencecprimeffectivemax')
       disp(shearconvergencecprimeffectivemax)
       disp('shearXcprimeffectivemax')
       disp(shearXcprimeffectivemax)
       disp('shearYcprimeffectivemax')
       disp(shearYcprimeffectivemax)
  
       
        disp('Hshearconvergencemax')
       disp(Hshearconvergencemax)
       disp('directionshearconvergencemax')
       disp(directionshearconvergencemax)
       disp('shearconvergenceshearconvergencemax')
       disp(shearconvergenceshearconvergencemax)
       disp('shearXshearconvergencemax')
       disp(shearXshearconvergencemax)
       disp('shearYshearconvergencemax')
       disp(shearYshearconvergencemax)
       
       disp('HconvergenceZmax')
       disp(HconvergenceZmax)
       disp('directionconvergenceZmax')
       disp(directionconvergenceZmax)
       disp('shearconvergenceconvergenceZmax')
       disp(shearconvergenceconvergenceZmax)
       disp('shearXconvergenceZmax')
       disp(shearXconvergenceZmax)
       disp('shearYconvergenceZmax')
       disp(shearYconvergenceZmax)
  
          disp('values at layer')
       
 disp('Cprimeffective')
       disp(Cprimeffective(645))
       disp('shearconvergence')
       disp(shearconvergence(645))
       disp('Convergencez')
       disp(convergenceZ(645))
       disp('shearx')
       disp(shearX(645))
       disp('sheary')
       disp(shearY(645))
%save HWM14our145500.mat VnX VnY VnZ ViZnop Uwind Ceffective Cprimeffective convergenceZ shearconvergence