clear all
close all
clc






%load MSISEDATA.mat
load constsandscaling.mat

%parameters for plotting

for h=1:Nz
    Dzzforplot(h)=Dzz(jmin,h)*(1000*1000/3600); %m^2/sec
    Dyyforplot(h)=Dyy(jmin,h)*(1000*1000/3600); %m^2/sec
    Dxxforplot(h)=Dxx(jmin,h)*(1000*1000/3600); %m^2/sec
    Dxzforplot(h)=Dxz(jmin,h)*(1000*1000/3600); %m^2/sec
    one(h)=1;
    niuinspmsiseforplot(h)= niuinspmsise(jmin,h)/3600; %sec^-1
    
    kappaspnormmsiseforplot(h)=kappaspnormmsise(jmin,h);
    
    VnXinViX(h)=(kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h)+COSI*COSI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
    VnYinViX(h)=(-kappaspnormmsiseforplot(h)*SINI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
    VnZinViX(h)=(-SINI*COSI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
    
    
    VnXinViY(h)=(kappaspnormmsiseforplot(h)*SINI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
    VnYinViY(h)=(kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h))/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
    VnZinViY(h)=(kappaspnormmsiseforplot(h)*COSI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));

    
    VnXinViZ(h)=(-SINI*COSI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
    VnYinViZ(h)=(-kappaspnormmsiseforplot(h)*COSI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
    VnZinViZ(h)=(kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h)+SINI*SINI)/(1+kappaspnormmsiseforplot(h)*kappaspnormmsiseforplot(h));
end
for h=2:Nz-1
    convergencecoeffx(h)=-((VnXinViZ(h+1)-VnXinViZ(h-1))/(z(h+1)-z(h-1)));
end
for h=2:Nz-1
    convergencecoeffy(h)=-((VnYinViZ(h+1)-VnYinViZ(h-1))/(z(h+1)-z(h-1)));
end

for h=2:Nz-1
    convergencecoeffx(h)=-((VnXinViZ(h+1)-VnXinViZ(h-1))/(z(h+1)-z(h-1)));
end

convergencecoeffx=smooth(convergencecoeffx);
convergencecoeffy=smooth(convergencecoeffy);

dphi=(0.25/9)*pi/10;
phi0=0;
Nphi=2*pi/dphi+1;
for p=1:Nphi
phi(p)=phi0+dphi*(p-1);
end
disp('phi')
disp(phi(215)*180/pi)
for h=2:Nz-1
    for p=1:Nphi
    Velcoef(h,p)=cos(phi(p))*VnXinViZ(h)+sin(phi(p))*VnYinViZ(h);
    convergencecoeff(h,p)=cos(phi(p))*convergencecoeffx(h)+sin(phi(p))*convergencecoeffy(h);
    
    end  
end
disp('cmax')
disp(max(max(Velcoef)))

for p=1:Nphi
    VelcoefPhi(p)=-Velcoef(600,p);
    convergencecoeffPhi(p)=convergencecoeff(600,p);
    nuli(p)=0;
end
for h=2:Nz-1
    
    Velcoef270(h)=-Velcoef(h,541);
    Velcoef275(h)=-Velcoef(h,551);
    Velcoef300(h)=-Velcoef(h,601);
    Velcoef330(h)=-Velcoef(h,661);
    Velcoef345(h)=-Velcoef(h,691);
    Velcoef360(h)=-Velcoef(h,721);
    
    
    
    convergencecoeff270(h)=convergencecoeff(h,541);
    convergencecoeff300(h)=convergencecoeff(h,601);
    convergencecoeff330(h)=convergencecoeff(h,691);
    convergencecoeff360(h)=convergencecoeff(h,721);
end
for h=1:Nz-1
    for p=1:Nphi
    Velcoef(h,p)=cos(phi(p))*VnXinViZ(h)+sin(phi(p))*VnYinViZ(h);
    %convergencecoeff(h,p)=cos(phi(p))*convergencecoeffx(h)+sin(phi(p))*convergencecoeffy(h);
    
    end  
end
disp(min(min(abs(Velcoef))))
for h=1:Nz-1
phiatbalancethor(h)=atan(-(VnXinViZ(h)/VnYinViZ(h)));
end

disp('phiatbalancethor')
disp(180+180*phiatbalancethor(735)/pi)

 for h=575:810
        phibalance11(h)=0;
          zbalance11(h)=0;
 end

 for p=1:round(0.5*Nphi)
 for h=575:810
    
        if   abs(Velcoef(h,p))<0.0001
     
          
          phibalance11(h)=phi(p);
          zbalance11(h)=z(h);
            
        end
end
 end
 k3=0;
 for h=575:810
     if phibalance11(h)>0 && zbalance11(h)>0
         k3=k3+1;
       phibalance1(k3)=phibalance11(h);
       
     zbalance1(k3)=zbalance11(h);
     end
 end
 
 for p=round(0.5*Nphi):Nphi
 for h=575:810
    
        if   abs(Velcoef(h,p))<0.0001
     
          
          phibalance11(h)=phi(p);
          zbalance11(h)=z(h);
            
        end
end
end
 k4=k3;
 for h=575:810
     if phibalance11(h)>0 && zbalance11(h)>0
         k4=k4+1;
       phibalance2(k4)=phibalance11(h);
       
     zbalance2(k4)=zbalance11(h);
     end
 end




p=1;
disp(phi(p)*(180/pi))
coefmax=max(convergencecoeff(2:Nz-1,p));
coefmin=min(convergencecoeff(2:Nz-1,p));

for k=2:Nz-1
   if convergencecoeff(k,p)==coefmax
       kmax0=k;
   end
end
disp('Kmax0=')
disp( kmax0)
disp('Zmax0=')
disp(Hnorm*z( kmax0))
disp('convergencecoeffmax0=')
disp(coefmax/(Hnorm*1000))

p=181;
disp(phi(p)*(180/pi))
coefmax=max(convergencecoeff(2:Nz-1,p));
coefmin=min(convergencecoeff(2:Nz-1,p));

for k=2:Nz-1
   if convergencecoeff(k,p)==coefmax
       kmax90=k;
   end
end
disp('Kmax90=')
disp( kmax90)
disp('Zmax90=')
disp(Hnorm*z( kmax90))
disp('convergencecoeffmax0=')
disp(coefmax/(Hnorm*1000))

p=181;
disp(phi(p)*(180/pi))
coefmax=max(convergencecoeff(2:Nz-1,p));
coefmin=min(convergencecoeff(2:Nz-1,p));

for k=2:Nz-1
   if convergencecoeff(k,p)==coefmin
       kmin90=k;
   end
end
disp('Kmin90=')
disp(kmin90)
disp('Zmin90=')
disp(Hnorm*z(kmin90))
disp('convergencecoeffmin90=')
disp(coefmin/(Hnorm*1000))
disp('velcoefmin90')
disp(Velcoef(kmin90,181))

p=361;
disp(phi(p)*(180/pi))
coefmax=max(convergencecoeff(2:Nz-1,p));
coefmin=min(convergencecoeff(2:Nz-1,p));

for k=2:Nz-1
   if convergencecoeff(k,p)==coefmin
       kmin180=k;
   end
end
disp('Kmin180=')
disp(kmin180)
disp('Zmin180=')
disp(Hnorm*z(kmin180))
disp('convergencecoeffmin180=')
disp(coefmin/(Hnorm*1000))
disp('velcoefmin180')
disp(Velcoef(kmin180,361))





p=541;
disp(phi(p)*(180/pi))
coefmax=max(convergencecoeff(2:Nz-1,p));
coefmin=min(convergencecoeff(2:Nz-1,p));

for k=2:Nz-1
   if convergencecoeff(k,p)==coefmax
       kmax270=k;
   end
end
disp('Kmax270=')
disp( kmax270)
disp('Zmax270=')
disp(Hnorm*z( kmax270))
disp('convergencecoeffmax0=')
disp(coefmax/(Hnorm*1000))

p=541;
disp(phi(p)*(180/pi))
coefmax=max(convergencecoeff(2:Nz-1,p));
coefmin=min(convergencecoeff(2:Nz-1,p));

for k=2:Nz-1
   if convergencecoeff(k,p)==coefmin
       kmin270=k;
   end
end
disp('Kmin270=')
disp(kmin270)
disp('Zmin270=')
disp(Hnorm*z(kmin270))
disp('convergencecoeffmin270=')
disp(coefmin/(Hnorm*1000))
disp('velcoefmin270')
disp(Velcoef(kmin270,181))

Nt=10500;
for tindex=1:round(Nt/50)
taunew(tindex)=tau(50*tindex);
Zmax0(tindex)= Hnorm*z(kmax0);
Zmax90(tindex)=Hnorm*z(kmax90);
Zmin90(tindex)=Hnorm*z(kmin90);
Zmin180(tindex)=Hnorm*z(kmin180);
Zmax270(tindex)=Hnorm*z(kmax270);
Zmin270(tindex)=Hnorm*z(kmin270);
end

kmin=375;
kmax=810;
for k=kmin+1:kmax-1
    for p=1:Nphi
 z1(k-kmin)=z(k);
 Velcoef1(k-kmin,p)=Velcoef(k,p);
 convergencecoeff1(k-kmin,p)=convergencecoeff(k,p);
    end
end

figure(1)

pcolor(phi(1:Nphi)*(180/pi), Hnorm*z1(2:kmax-kmin-1),  -Velcoef1(2:kmax-kmin-1,1:Nphi));
hold on
%grid(gca,'minor')
%view(90,-90);
 title ('Ion vertical drift factor','FontSize',16,'FontWeight','bold')
 set(gca,'FontSize',16,'Fontweight','bold')
 ylabel('Height h, km','Fontsize',16,'Fontweight','bold');
 xlabel('Wind direction angle, degrees','Fontsize',16,'Fontweight','bold');
 xticks([0 45 90 135 180 225 270 315 360])
 axis([-3,360,74,Hnorm*z1(kmax-kmin-1)])
 yticks([80 100 120 140])
 
 %%%Long arrows left
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [195-180 195-180];          % the location in data units
ha.Y = [140 120];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;



ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [235-180 235-180];          % the location in data units
ha.Y = [140 125]; 
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [255-180 255-180];          % the location in data units
ha.Y = [140 130];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275-180 275-180];          % the location in data units
ha.Y = [140 135];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%Small arrows at balance1
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335-180 335-180];          % the location in data units
ha.Y = [119 124];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335-180 335-180];          % the location in data units
ha.Y = [117 112];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350-180 350-180];          % the location in data units
ha.Y = [113 118];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350-180 350-180];          % the location in data units
ha.Y = [111 106];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;




ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320-180 320-180];          % the location in data units
ha.Y = [123 128];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320-180 320-180];          % the location in data units
ha.Y = [121 116];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307-180 307-180];          % the location in data units
ha.Y = [126 121];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307-180 307-180];          % the location in data units
ha.Y = [128 133];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295-180 295-180];          % the location in data units
ha.Y = [133 138];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295-180 295-180];          % the location in data units
ha.Y = [131 126];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


%Long arrows middle
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [115 115];          % the location in data units
ha.Y = [140 145];

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;
 
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [135 135];          % the location in data units
ha.Y = [140 150];
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;
 

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [155 155];          % the location in data units
ha.Y = [140 155];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;
 

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [195 195];          % the location in data units
ha.Y = [140 160];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [235 235];          % the location in data units
ha.Y = [140 155];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [255 255];          % the location in data units
ha.Y = [140 150];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275 275];          % the location in data units
ha.Y = [140 145];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


%Small arrows for balance 2
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335 335];          % the location in data units
ha.Y = [124 119];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;




ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335 335];          % the location in data units
ha.Y = [112 117];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350 350];          % the location in data units
ha.Y = [118 113];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350 350];          % the location in data units
ha.Y = [106 111];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [128 123];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [116 121];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307 307];          % the location in data units
ha.Y = [121 126];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307 307];          % the location in data units
ha.Y = [133 128];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295 295];          % the location in data units
ha.Y = [138 133];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295 295];          % the location in data units
ha.Y = [126 131];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%%Long arrows right
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [140 132];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335 335];          % the location in data units
ha.Y = [140 127];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350 350];          % the location in data units
ha.Y = [140 122];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%Arrows for contour 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [185-180 185-180];          % the location in data units
ha.Y = [108.2000 108.2000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [200-180 200-180];          % the location in data units
ha.Y = [104 104-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [215-180 215-180];          % the location in data units
ha.Y = [101.8000 101.8000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [230-180 230-180];          % the location in data units
ha.Y = [100.4000 100.4000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [245-180 245-180];          % the location in data units
ha.Y = [99.4000 99.4000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [260-180 260-180];          % the location in data units
ha.Y = [99 99-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275-180 275-180];          % the location in data units
ha.Y = [99 99-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [290-180 290-180];          % the location in data units
ha.Y = [99.4000 99.4000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [305-180 305-180];          % the location in data units
ha.Y = [100.2000 100.2000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320-180 320-180];          % the location in data units
ha.Y = [101.6000 101.6000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [334-180 334-180];          % the location in data units
ha.Y = [103.8000 103.8000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%Arrows for contour 2

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [185 185];          % the location in data units
ha.Y = [108.2000 108.2000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [200 200];          % the location in data units
ha.Y = [104 104+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [215 215];          % the location in data units
ha.Y = [101.8000 101.8000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [230 230];          % the location in data units
ha.Y = [100.4000 100.4000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [245 245];          % the location in data units
ha.Y = [99.4000 99.4000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [260 260];          % the location in data units
ha.Y = [99 99+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275 275];          % the location in data units
ha.Y = [99 99+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [290 290];          % the location in data units
ha.Y = [99.4000 99.4000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [305 305];          % the location in data units
ha.Y = [100.2000 100.2000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [101.6000 101.6000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [334 334];          % the location in data units
ha.Y = [103.8000 103.8000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('line');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [98 103];          % the location in data units
ha.Y = [155 155];  
ha.LineWidth  = 1;
ha = annotation('line');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [98+180 103+180];          % the location in data units
ha.Y = [155 155];  
ha.LineWidth  = 1;
 colorbar;
%colorbar;
 fclose all;
shading interp;
%grid on
hold on
contour(phi(1:Nphi)*(180/pi), Hnorm*z1(2:kmax-kmin-1),  -Velcoef1(2:kmax-kmin-1,1:Nphi))
hold on
plot(180+180*phiatbalancethor(500:800)/pi,z(500:800)*Hnorm,'b','linewidth',1)
hold on
plot(360+180*phiatbalancethor(500:800)/pi,z(500:800)*Hnorm,'y','linewidth',1)
hold on
text(40,80, '(a)', 'Fontsize',14, 'color', 'w', 'fontweight','bold')
text(103,155, 'C=0', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
text(103+180,155, 'C=0', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
hold off

set(figure(1), 'PaperPositionMode', 'manual');
set(figure(1), 'PaperUnits', 'centimeters');
set(figure(1), 'Paperposition', [0 0 18 8]);


print(figure(1),'Fig1-1','-djpeg', '-r100')







figure(2)

pcolor(phi(1:Nphi)*(180/pi), Hnorm*z1(2:kmax-kmin-1),  convergencecoeff1(2:kmax-kmin-1,1:Nphi)/(20000));
hold on
%grid(gca,'minor')
%view(90,-90);
 title ('Ion convergence/divergence factor, m^{-1}','FontSize',16,'FontWeight','bold')
 set(gca,'FontSize',16,'Fontweight','bold')
 ylabel('Height h, km','Fontsize',16,'Fontweight','bold');
 xlabel('Wind direction angle, degrees','Fontsize',16,'Fontweight','bold');
 xticks([0 45 90 135 180 225 270 315 360])
 axis([-3,360,74,Hnorm*z1(kmax-kmin-1)])
 yticks([80 100 120 140])
 
 %%%Long arrows left
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [195-180 195-180];          % the location in data units
ha.Y = [140 120];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;



ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [235-180 235-180];          % the location in data units
ha.Y = [140 125]; 
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [255-180 255-180];          % the location in data units
ha.Y = [140 130];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275-180 275-180];          % the location in data units
ha.Y = [140 135];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%Small arrows at balance1
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335-180 335-180];          % the location in data units
ha.Y = [119 124];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335-180 335-180];          % the location in data units
ha.Y = [117 112];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350-180 350-180];          % the location in data units
ha.Y = [113 118];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350-180 350-180];          % the location in data units
ha.Y = [111 106];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;




ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320-180 320-180];          % the location in data units
ha.Y = [123 128];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320-180 320-180];          % the location in data units
ha.Y = [121 116];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307-180 307-180];          % the location in data units
ha.Y = [126 121];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307-180 307-180];          % the location in data units
ha.Y = [128 133];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295-180 295-180];          % the location in data units
ha.Y = [133 138];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295-180 295-180];          % the location in data units
ha.Y = [131 126];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


%Long arrows middle
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [115 115];          % the location in data units
ha.Y = [140 145];

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;
 
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [135 135];          % the location in data units
ha.Y = [140 150];
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;
 

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [155 155];          % the location in data units
ha.Y = [140 155];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;
 

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [195 195];          % the location in data units
ha.Y = [140 160];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [235 235];          % the location in data units
ha.Y = [140 155];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [255 255];          % the location in data units
ha.Y = [140 150];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275 275];          % the location in data units
ha.Y = [140 145];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


%Small arrows for balance 2
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335 335];          % the location in data units
ha.Y = [124 119];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;




ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335 335];          % the location in data units
ha.Y = [112 117];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350 350];          % the location in data units
ha.Y = [118 113];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350 350];          % the location in data units
ha.Y = [106 111];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [128 123];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [116 121];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307 307];          % the location in data units
ha.Y = [121 126];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [307 307];          % the location in data units
ha.Y = [133 128];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295 295];          % the location in data units
ha.Y = [138 133];   

ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [295 295];          % the location in data units
ha.Y = [126 131];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%%Long arrows right
ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [140 132];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [335 335];          % the location in data units
ha.Y = [140 127];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [350 350];          % the location in data units
ha.Y = [140 122];   
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%Arrows for contour 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [185-180 185-180];          % the location in data units
ha.Y = [108.2000 108.2000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [200-180 200-180];          % the location in data units
ha.Y = [104 104-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [215-180 215-180];          % the location in data units
ha.Y = [101.8000 101.8000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [230-180 230-180];          % the location in data units
ha.Y = [100.4000 100.4000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [245-180 245-180];          % the location in data units
ha.Y = [99.4000 99.4000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [260-180 260-180];          % the location in data units
ha.Y = [99 99-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275-180 275-180];          % the location in data units
ha.Y = [99 99-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [290-180 290-180];          % the location in data units
ha.Y = [99.4000 99.4000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [305-180 305-180];          % the location in data units
ha.Y = [100.2000 100.2000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320-180 320-180];          % the location in data units
ha.Y = [101.6000 101.6000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [334-180 334-180];          % the location in data units
ha.Y = [103.8000 103.8000-5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

%Arrows for contour 2

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [185 185];          % the location in data units
ha.Y = [108.2000 108.2000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [200 200];          % the location in data units
ha.Y = [104 104+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [215 215];          % the location in data units
ha.Y = [101.8000 101.8000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [230 230];          % the location in data units
ha.Y = [100.4000 100.4000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [245 245];          % the location in data units
ha.Y = [99.4000 99.4000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [260 260];          % the location in data units
ha.Y = [99 99+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [275 275];          % the location in data units
ha.Y = [99 99+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [290 290];          % the location in data units
ha.Y = [99.4000 99.4000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [305 305];          % the location in data units
ha.Y = [100.2000 100.2000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;

ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [320 320];          % the location in data units
ha.Y = [101.6000 101.6000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('arrow');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [334 334];          % the location in data units
ha.Y = [103.8000 103.8000+5];  
ha.LineWidth  = 1;          % make the arrow bolder for the picture
ha.HeadWidth  = 5;
ha.HeadLength = 5;


ha = annotation('line');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [98 103];          % the location in data units
ha.Y = [155 155];  
ha.LineWidth  = 1;
ha = annotation('line');  % store the arrow information in ha
ha.Parent = gca;           % associate the arrow the the current axes
ha.X = [98+180 103+180];          % the location in data units
ha.Y = [155 155];  
ha.LineWidth  = 1;
 colorbar;
%colorbar;
 fclose all;
shading interp;
%grid on
hold on
contour(phi(1:Nphi)*(180/pi), Hnorm*z1(2:kmax-kmin-1),  convergencecoeff1(2:kmax-kmin-1,1:Nphi)/(20000))
hold on
plot(180+180*phiatbalancethor(500:800)/pi,z(500:800)*Hnorm,'b','linewidth',1)
hold on
plot(360+180*phiatbalancethor(500:800)/pi,z(500:800)*Hnorm,'y','linewidth',1)
hold on
text(40,80, '(b)', 'Fontsize',14, 'color', 'w', 'fontweight','bold')
text(103,155, 'C=0', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
text(103+180,155, 'C=0', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
hold off

set(figure(2), 'PaperPositionMode', 'manual');
set(figure(2), 'PaperUnits', 'centimeters');
set(figure(2), 'Paperposition', [0 0 18 8]);


print(figure(2),'Fig1-2','-djpeg', '-r100')




%
%{
Nt=3250;

hold on


figure(3)
%subplot(1,3,1)


load ZTU0=0-Da-z095forfigsj0281.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 
  set(gca,'FontSize',16,'Fontweight','bold')
  
 % title ('Electron density Ne(h,t)/Nom','Fontsize',10,'FontWeight','bold')
% xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 colorbar;
 %cbh=colorbar('v');
 %caxis([0 1])
%set(cbh,'XTick',[0:0.2:1])
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, 'A1', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
 axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off
set(figure(3), 'PaperPositionMode', 'manual');
set(figure(3), 'PaperUnits', 'centimeters');
set(figure(3), 'Paperposition', [0 0 7 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(3),'Fig2-1','-dpng', '-r100','-r100' )


%subplot(1,3,2)
figure(4)
load ZTU0=0-Da-z0120forfigsj0281.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);

  set(gca,'FontSize',16,'Fontweight','bold')
  
% xlabel('Time, hours');
 %ylabel('Height h, km');
%cbh=colorbar('v');
 %caxis([0 1])
%set(cbh,'XTick',[0:0.2:1])
colorbar
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, 'B1', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
 
hold off
set(figure(4), 'PaperPositionMode', 'manual');
set(figure(4), 'PaperUnits', 'centimeters');
set(figure(4), 'Paperposition', [0 0 6 7]);
%set(colorbar,'position', [0.9 0.2 0.05 0.7])
%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(4),'Fig2-2','-dpng', '-r100','-r100')

%subplot(1,3,3)
figure(5)


load ZTU0=0-Da-z0135forfigsj0281.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 
  set(gca,'FontSize',16,'Fontweight','bold')
  %title ('Electron density Ne(h,t)/Nom','Fontsize',10,'FontWeight','bold')
 %xlabel('Time, hours');
 %ylabel('Height h, km');
 %cbh=colorbar('v');
 %caxis([0 1])
%set(cbh,'XTick',[0:0.2:1])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, 'C1', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
%text(-1, 164,'Electron density N_{e}(h,t)/N_{om}','Fontsize',10,'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

set(figure(5), 'PaperPositionMode', 'manual');
set(figure(5), 'PaperUnits', 'centimeters');
set(figure(5), 'Paperposition', [0 0 6 7]);



print(figure(5),'Fig2-3','-dpng', '-r100')
figure(6)

%subplot(1,3,1)


load ZTU0=0-Da-z0=95-Heff=8j0281.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
cbh=colorbar('v');
 caxis([0 1])
set(cbh,'XTick',[0:0.2:1])
%colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, 'A2', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off
%text(65,0.2,'Time, hours','FontSize',16,'Fontweight','bold')
set(figure(6), 'PaperPositionMode', 'manual');
set(figure(6), 'PaperUnits', 'centimeters');
set(figure(6), 'Paperposition', [0 0 7 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(6),'Fig2-4','-dpng', '-r100')

%subplot(1,3,2)
figure(7)
load ZTU0=0-Da-z0=120-Heff=8j0281.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 %ylabel('Height h, km');
 cbh=colorbar('v');
 caxis([0 1])
set(cbh,'XTick',[0:0.2:1])
%colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, 'B2', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
 %text(65,0.2,'Time, hours','FontSize',16,'Fontweight','bold')
hold off

set(figure(7), 'PaperPositionMode', 'manual');
set(figure(7), 'PaperUnits', 'centimeters');
set(figure(7), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(7),'Fig2-5','-dpng', '-r100')

%subplot(1,3,3)
figure(8)

load ZTU0=0-Da-z0=135-Heff=8j0281.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
% xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 %ylabel('Height h, km');
 cbh=colorbar('v');
 caxis([0 1])
set(cbh,'XTick',[0:0.2:1])
%colorbar
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, 'C2', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
 %text(65,0.2,'Time, hours','FontSize',16,'Fontweight','bold')
hold off

set(figure(8), 'PaperPositionMode', 'manual');
set(figure(8), 'PaperUnits', 'centimeters');
set(figure(8), 'Paperposition', [0 0 6 7]);


print(figure(8),'Fig2-6','-dpng', '-r100')



%{
subplot(1,4,4)
% plot(z*Hnorm, Dxxforplot)
 
%hold on
 %plot(z*Hnorm, Dyyforplot)

 plot(z*Hnorm, Dzzforplot)
 view(90,-90)
 hold on
 %plot(z*Hnorm, Dxzforplot)
 xlabel('Height, km');
 %ylabel('Diffusion coefficient, m^2 sec^-1');
 %legend('D_{zz}')
 %xlabel('Height, km');
 title('Diffusion coefficient, m^2 sec^-1');
%axis([90,200,0,0.5*10^4])
yticks([0 0.1 1 10 100 1000 10000 100000])
 %}


figure(9)



%subplot(1,3,1)


load ZTU0=25-Da-alpha0-z095forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 
  set(gca,'FontSize',16,'Fontweight','bold')
  %title ('Electron density Ne(h,t)/Nom','FontSize',8,'FontWeight','bold')
 %xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 %cbh=colorbar('v');
  %caxis([0 2])
%set(cbh,'XTick',[0:0.2:1])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax0,'y','linewidth',1)
text(0.03,155, 'A1-N', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off
set(figure(9), 'PaperPositionMode', 'manual');
set(figure(9), 'PaperUnits', 'centimeters');
set(figure(9), 'Paperposition', [0 0 7 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(9),'Fig3-1','-dpng', '-r100')


%subplot(1,3,2)
figure(10)
Nt=Nt;
for tindex=1:round(Nt/50)
taunew(tindex)=tau(50*tindex);
end

load ZTU0=25-Da-alpha0-z0120forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);

  set(gca,'FontSize',16,'Fontweight','bold')

  % title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
   
% xlabel('Time, hours');
 %ylabel('Height h, km');
  %caxis([0 2])
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.4:2])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax0,'y','linewidth',1)
text(0.03,155, 'B1-N', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off
set(figure(10), 'PaperPositionMode', 'manual');
set(figure(10), 'PaperUnits', 'centimeters');
set(figure(10), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(10),'Fig3-2','-dpng', '-r100')


%subplot(1,3,3)
figure(11)
Nt=Nt;
for tindex=1:round(Nt/50)
taunew(tindex)=tau(50*tindex);
end

load ZTU0=25-Da-alpha0-z0135forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 
  set(gca,'FontSize',16,'Fontweight','bold')
  %title ('Electron density Ne(h,t)/Nom','FontSize',8,'FontWeight','bold')
 %xlabel('Time, hours');
 %ylabel('Height h, km');
  %caxis([0 2])

 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.4:2])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax0,'y','linewidth',1)
text(0.03,155, 'C1-N', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

set(figure(11), 'PaperPositionMode', 'manual');
set(figure(11), 'PaperUnits', 'centimeters');
set(figure(11), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(11),'Fig3-3','-dpng', '-r100')

figure(12)
%subplot(1,3,1)


load ZTU0=25-Da-alpha180-z095forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
  %caxis([0 1.2])
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1.2])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmin180,'blue','linewidth',1)
text(0.03,155, 'A1-S', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

set(figure(12), 'PaperPositionMode', 'manual');
set(figure(12), 'PaperUnits', 'centimeters');
set(figure(12), 'Paperposition', [0 0 7 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(12),'Fig3-4','-dpng', '-r100')

%subplot(1,3,2)

figure(13)


load ZTU0=25-Da-alpha180-z0120forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 %ylabel('Height h, km');
  %caxis([0 1.2])
%cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1.2])
colorbar
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmin180,'blue','linewidth',1)
text(0.03,155, 'B1-S', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

set(figure(13), 'PaperPositionMode', 'manual');
set(figure(13), 'PaperUnits', 'centimeters');
set(figure(13), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(13),'Fig3-5','-dpng', '-r100')
%subplot(1,3,3)

figure(14)
load ZTU0=25-Da-alpha180-z0135forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 %ylabel('Height h, km');
  %caxis([0 1.2])
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1.2])
colorbar
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmin180,'blue','linewidth',1)
text(0.03,155, 'C1-S', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

 %suptitle('Electron density N_{e}(h,t)/N_{om}')
set(figure(14), 'PaperPositionMode', 'manual');
set(figure(14), 'PaperUnits', 'centimeters');
set(figure(14), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(14),'Fig3-6','-dpng', '-r100')
%{
subplot(1,4,4)
plot(z*Hnorm, -VnXinViZ,'-b')
hold on
plot(z(4:1:Nz-3)*Hnorm, smooth(convergencecoeffx(4:1:Nz-3)),':g')
hold on
%plot(z*Hnorm, VnZinViZ,'-g')
view(90,-90)
xlabel('Height, km');
title('Velocity and convergence factors')
%legend('Factor of VnX in ViZ', 'Factor of VnY in ViZ', 'Factor of VnZ in ViZ' )
%}

figure(15)

%subplot(1,3,1)


load ZTU0=25-Da-alpha90-z095forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
  %caxis([0 2])
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax90,'y',taunew, Zmin90,'blue','linewidth',1)
text(0.03,155, 'A1-W', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off
%suptitle('Electron density N_{e}(h,t)/N_{om}')
set(figure(15), 'PaperPositionMode', 'manual');
set(figure(15), 'PaperUnits', 'centimeters');
set(figure(15), 'Paperposition', [0 0 7 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(15),'Fig4-1','-dpng', '-r100')

%subplot(1,3,2)

figure(16)
load ZTU0=25-Da-alpha90-z0120forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 
  set(gca,'FontSize',16,'Fontweight','bold')
  %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
 %xlabel('Time, hours');
 %ylabel('Height h, km');
 % caxis([0 2])
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.4:2])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax90,'y',taunew, Zmin90,'blue','linewidth',1)
text(0.03,155, 'B1-W', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off
set(figure(16), 'PaperPositionMode', 'manual');
set(figure(16), 'PaperUnits', 'centimeters');
set(figure(16), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(16),'Fig4-2','-dpng', '-r100')
%subplot(1,3,3)
figure(17)


load ZTU0=25-Da-alpha90-z0135forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours');
 %ylabel('Height h, km');
   %caxis([0 1])
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax90,'y',taunew, Zmin90,'blue','linewidth',1)
text(0.03,155, 'C1-W', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

set(figure(17), 'PaperPositionMode', 'manual');
set(figure(17), 'PaperUnits', 'centimeters');
set(figure(17), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(17),'Fig4-3','-dpng', '-r100')

figure(18)
%subplot(1,3,1)


load ZTU0=25-Da-alpha270-z095forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
% xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 %caxis([0 1.3])
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1.5])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax270,'y',taunew, Zmin270,'blue','linewidth',1)
text(0.03,155, 'A1-E', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off
set(figure(18), 'PaperPositionMode', 'manual');
set(figure(18), 'PaperUnits', 'centimeters');
set(figure(18), 'Paperposition', [0 0 7 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(18),'Fig4-4','-dpng', '-r100')

%subplot(1,3,2)

figure(19)
load ZTU0=25-Da-alpha270-z0120forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 %ylabel('Height h, km');
%cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1.2])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax270,'y',taunew, Zmin270,'blue','linewidth',1)
text(0.03,155, 'B1-E', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

set(figure(19), 'PaperPositionMode', 'manual');
set(figure(19), 'PaperUnits', 'centimeters');
set(figure(19), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(19),'Fig4-5','-dpng', '-r100')

%subplot(1,3,3)
figure(20)

load ZTU0=25-Da-alpha270-z0135forfigj0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
% xlabel('Time, hours','FontSize',16,'Fontweight','bold');
% ylabel('Height h, km');
%cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:1.2])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
hold on
plot(taunew, Zmax270,'y',taunew, Zmin270,'blue','linewidth',1)
text(0.03,155, 'C1-E', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,0.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=[0 0.2 0.5];
 xticks(timeticks)
hold off

set(figure(20), 'PaperPositionMode', 'manual');
set(figure(20), 'PaperUnits', 'centimeters');
set(figure(20), 'Paperposition', [0 0 6 7]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(20),'Fig4-6','-dpng', '-r100')
%{
subplot(1,4,4)

%hold on
plot(z(4:1:Nz-3)*Hnorm, smooth(convergencecoeffy(4:1:Nz-3)),':g')
hold on
plot(z*Hnorm, -VnYinViZ,'-b')
hold on
view(90,-90)
xlabel('Height, km');
title('Velocity and convergence factors')
%legend('Factor of VnX in ViZ', 'Factor of VnY in ViZ', 'Factor of VnZ in ViZ' )

%}


%{
subplot(1,4,4)
%hold on
plot(z(4:1:Nz-3)*Hnorm, -smooth(convergencecoeffx(4:1:Nz-3)),':g')
hold on
plot(z*Hnorm, VnXinViZ,'-b')
hold on
view(90,-90)
xlabel('Height, km');
title('Velocity and convergence factors')
%}


%{
subplot(1,4,4)
%hold on
plot(z(4:1:Nz-3)*Hnorm, -smooth(convergencecoeffy(4:1:Nz-3)),':g')
hold on
plot(z*Hnorm, VnYinViZ,'-b')
hold on
view(90,-90)
xlabel('Height, km');
title('Velocity and convergence factors')
%legend('Factor of VnX in ViZ', 'Factor of VnY in ViZ', 'Factor of VnZ in ViZ' )
%}

 

%{
figure(8)

plot(convergencecoeff330(2:Nz-1),Hnorm*z(2:Nz-1),'b',convergencecoeff330(2:Nz-1),Hnorm*z(2:Nz-1),'g')

figure(9)
plot(z*Hnorm, kappaspnormmsiseforplot)
hold on
view(90,-90)
plot(z*Hnorm, one,'--')
 xlabel('Height, km','FontSize',14);
 ylabel('k=\nu_{in}/\omega','FontSize',14);
 axis([40,200,0,10])
 
 
 figure(10)
  plot(z*Hnorm, Dxxforplot)
 
hold on
 plot(z*Hnorm, Dyyforplot)

 plot(z*Hnorm, Dzzforplot)
 view(90,-90)
 hold on
 plot(z*Hnorm, Dxzforplot)
 xlabel('Height, km');
 ylabel('Diffusion coefficient, m^2 sec^-1');
 legend('D_{zz}')
 xlabel('Height, km');
 title('Diffusion coefficient, m^2 sec^-1');
axis([90,200,0,0.5*10^4])
yticks([0 0.1 1 10 100 1000 10000 100000])
 
dHeff=0.2;
Heffmin=8;
Heffmax=35;
NHeff=(Heffmax-Heffmin)/dHeff+1;

p=1;
zom=z(475);
for i=1:NHeff
Heff(i)=(Heffmin+dHeff*(i-1))/Hnorm;
end
for k=2:Nz-2
   for i=1:NHeff
    Concentrationrate(k,i)=convergencecoeff(k,p)+2*Velcoef(k,p)*((z(k)-zom)/(Heff(i)*Heff(i)));
   end
end
for k=2:Nz-2
    Concentrationrate1D(k)=Concentrationrate(k,20);
end
figure(11)

pcolor(Hnorm*Heff(1:NHeff), Hnorm*z(kmin+1:kmax-2) , Concentrationrate(kmin+1:kmax-2,1:NHeff));
hold on
%grid(gca,'minor')
%view(90,-90);
 title ('Concentrationrate','FontSize',16,'FontWeight','bold')
 ylabel('Height, km','FontSize',14);
 xlabel('Heff','FontSize',14);
  colorbar;
%colorbar;
 fclose all;
shading interp;
%grid off
hold on
contour(Hnorm*Heff(1:NHeff), Hnorm*z(kmin+1:kmax-2) , Concentrationrate(kmin+1:kmax-2,1:NHeff))
hold off

figure(12)
view(-90,90)
hold on
plot( Hnorm*z(2:Nz-2), Concentrationrate1D(2:Nz-2))
%}
Nt=6500;
for tindex=1:round(Nt/50)
taunew(tindex)=tau(50*tindex);
end
figure(21)

%subplot(1,2,1)

%}

for tindex=1:round(Nt/50)
taunew(tindex)=tau(50*tindex);
end
disp(taunew(Nt/50))
load ZTU0=150-Da-alpha0-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
figure(21)
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);

  set(gca,'FontSize',16,'Fontweight','bold')
  
% title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold','position', [1.7 ,160])
 %xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
%cbh=colorbar('v');
%set(cbh,'XTick',[0:1:6])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(n)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
 yticks([80 100 120 140])
hold off
set(figure(21), 'PaperPositionMode', 'manual');
set(figure(21), 'PaperUnits', 'centimeters');
set(figure(21), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(21),'Fig5 for 150-1','-dpng', '-r100')

%subplot(1,2,2)
figure(22)
load ZTU0=150-Da-alpha45-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours');
 %ylabel('Height h, km');
cbh=colorbar('v');
set(cbh,'XTick',[0.2 2 4 6])
%colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(n-w)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off
set(figure(22), 'PaperPositionMode', 'manual');
set(figure(22), 'PaperUnits', 'centimeters');
set(figure(22), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(22),'Fig5 for 150-2','-dpng', '-r100')

figure(23)

%subplot(1,2,1)


load ZTU0=150-Da-alpha90-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 cbh=colorbar('v');
set(cbh,'XTick',[0.2 2 4 6])
 
 %colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(w)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off
set(figure(23), 'PaperPositionMode', 'manual');
set(figure(23), 'PaperUnits', 'centimeters');
set(figure(23), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(23),'Fig5 for 150-3','-dpng', '-r100')

%subplot(1,2,2)
figure(24)

load ZTU0=150-Da-alpha135-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
% xlabel('Time, hours');
 %ylabel('Height h, km');
 colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(w-s)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off

set(figure(24), 'PaperPositionMode', 'manual');
set(figure(24), 'PaperUnits', 'centimeters');
set(figure(24), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(24),'Fig5 for 150-4','-dpng', '-r100')


figure(25)

%subplot(1,2,1)


load ZTU0=150-Da-alpha180-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
% xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:9])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(s)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off
set(figure(25), 'PaperPositionMode', 'manual');
set(figure(25), 'PaperUnits', 'centimeters');
set(figure(25), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(25),'Fig5 for 150-5','-dpng', '-r100')

%subplot(1,2,2)

figure(26)
load ZTU0=150-Da-alpha225-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours');
 %ylabel('Height h, km');
 %cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:0.9])
colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(s-e)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off

set(figure(26), 'PaperPositionMode', 'manual');
set(figure(26), 'PaperUnits', 'centimeters');
set(figure(26), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(26),'Fig5 for 150-6','-dpng', '-r100')





figure(27)
load ZTU0=150-Da-alpha270-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours');
 ylabel('Height h, km');
 cbh=colorbar('v');
%set(cbh,'XTick',[0:0.2:0.9])
%cbh=colorbar('v');
 %caxis([0 1])
set(cbh,'XTick',[0:0.3:1.2])
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(e)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off

set(figure(27), 'PaperPositionMode', 'manual');
set(figure(27), 'PaperUnits', 'centimeters');
set(figure(27), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(27),'Fig5 for 150-7','-dpng', '-r100')


%subplot(1,2,2)

figure(28)
load ZTU0=150-Da-alpha315-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 %ylabel('Height h, km');
 colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(e-n)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])

 
 hold off
set(figure(28), 'PaperPositionMode', 'manual');
set(figure(28), 'PaperUnits', 'centimeters');
set(figure(28), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(28),'Fig5 for 150-8','-dpng', '-r100')


figure(29)

%subplot(1,2,1)


load ZTU0=100-Da-alpha0-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);

  set(gca,'FontSize',16,'Fontweight','bold')
   %title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold','position', [1.7 ,160])
 %xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(n)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off
set(figure(29), 'PaperPositionMode', 'manual');
set(figure(29), 'PaperUnits', 'centimeters');
set(figure(29), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(29),'Fig6-1','-dpng', '-r100')
%subplot(1,2,2)
figure(30)
load ZTU0=100-Da-alpha45-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;

kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Electron density Ne(h,t)/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')

  %xlabel('Time, hours');
 %ylabel('Height h, km');
 cbh=colorbar('v');
 %caxis([0 10])
 set(cbh,'XTick',[0.2 2 4 6])
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(n-w)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off
 
set(figure(30), 'PaperPositionMode', 'manual');
set(figure(30), 'PaperUnits', 'centimeters');
set(figure(30), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(30),'Fig6-2','-dpng', '-r100')

figure(31)

%subplot(1,2,1)


load ZTU0=100-Da-alpha90-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
% xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
  cbh=colorbar('v');
 %caxis([0 10])
 set(cbh,'XTick',[0.2 2 4 6])
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(w)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off

set(figure(31), 'PaperPositionMode', 'manual');
set(figure(31), 'PaperUnits', 'centimeters');
set(figure(31), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(31),'Fig6-3','-dpng', '-r100')
%subplot(1,2,2)
figure(32)

load ZTU0=100-Da-alpha135-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
% title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')

  %xlabel('Time, hours');
 %ylabel('Height h, km');
 cbh=colorbar('v');
set(cbh,'XTick',[0.1 1 2 3])
 %colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(w-s)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off


set(figure(32), 'PaperPositionMode', 'manual');
set(figure(32), 'PaperUnits', 'centimeters');
set(figure(32), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(32),'Fig6-4','-dpng', '-r100')
figure(33)

%subplot(1,2,1)


load ZTU0=100-Da-alpha180-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 
  %xlabel('Time, hours');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(s)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off

set(figure(33), 'PaperPositionMode', 'manual');
set(figure(33), 'PaperUnits', 'centimeters');
set(figure(33), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(33),'Fig6-5','-dpng', '-r100')
%subplot(1,2,2)

figure(34)
load ZTU0=100-Da-alpha225-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
 set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',14);
 %ylabel('Height h, km','FontSize',14);
 colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(s-e)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off


set(figure(34), 'PaperPositionMode', 'manual');
set(figure(34), 'PaperUnits', 'centimeters');
set(figure(34), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(34),'Fig6-6','-dpng', '-r100')
figure(35)

%subplot(1,2,1)


load ZTU0=100-Da-alpha270-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 ylabel('Height h, km','FontSize',16,'Fontweight','bold');
 
 cbh=colorbar('v');
 %caxis([0 10])
 set(cbh,'XTick',[0.3:0.3:1.2])
 %colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(e)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off

set(figure(35), 'PaperPositionMode', 'manual');
set(figure(35), 'PaperUnits', 'centimeters');
set(figure(35), 'Paperposition', [0 0 10 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(35),'Fig6-7','-dpng', '-r100')

%subplot(1,2,2)

figure(36)
load ZTU0=100-Da-alpha315-z0=120-Heff=8j0221.mat
%display (size(NZT(5,:)))
t=50:50:Nt;
kmin=375;
kmax=800;
pcolor( taunew(1:round(Nt/50)), Hnorm*z1(1:kmax-kmin-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)));
hold on
%grid(gca,'minor')
%view(0,-90);
 %title ('Normalized electron density Ne/Nom','FontSize',16,'FontWeight','bold')
  set(gca,'FontSize',16,'Fontweight','bold')
 %xlabel('Time, hours','FontSize',16,'Fontweight','bold');
 %ylabel('Height h, km');
 colorbar;
 fclose all;
shading interp;
%grid off
%colorbar;
hold on
contour(taunew(1:round(Nt/50)), Hnorm*z(kmin+1:kmax-1), NZT(1+175:kmax-kmin-1+175,1:round(Nt/50)))
text(0.03,155, '(e-n)', 'Fontsize',12, 'color', 'w', 'fontweight','bold')
axis([-0.005,1.5,74,Hnorm*z1(kmax-kmin-1)])
 timeticks=0:0.4:1.4;
 xticks(timeticks)
  yticks([80 100 120 140])
hold off

set(figure(36), 'PaperPositionMode', 'manual');
set(figure(36), 'PaperUnits', 'centimeters');
set(figure(36), 'Paperposition', [0 0 9 5]);

%suptitle('Electron density N_{e}(h,t)/N_{om}')
print(figure(36),'Fig6-8','-dpng', '-r100')


%}
disp(taunew(round(Nt/50)))