% Plots density contours from the spherical grid

clf

% Read the density

path0 = './IonEnergyDeposition.dat'
path1 = './IonEnergyDeposition-1s.dat'
path2 = './IonEnergyDeposition-2s.dat'
path3 = './IonEnergyDeposition-3s.dat'


grdden = load(path0);
grdden1 = load(path1);
grdden2 = load(path2);
grdden3 = load(path3);



kedep = grdden(:,1);   %+grdden2(:,1)+grdden3(:,1);
%kedep = grdden1(:,1)+grdden2(:,1)+grdden3(:,1);

% number of cells

  nr  = 71;
  nth = 36;
  nph = 73;

  
% Cell size

pi = acos(-1.0)

rlen = 450. - 100.;
thlen = pi - 5.0*pi/180.;
thb1 = 2.5*pi/180.;
thb2 = pi - thb1;
phlen = 2.*acos(-1.0);

dr = rlen/(nr-1);
dth = thlen/(nth-1);
dph = phlen/(nph-1);

[th1,ph1,r1] = meshgrid(thb1:dth:(nth)*dth, 0:dph:(nph-1)*dph,...
                 0:dr:(nr-1)*dr);

%[ph1,th1,r1] = ndgrid(0:dph:(nph-1)*dph, thb1:dth:(nth)*dth, 0:dr:(nr-1)*dr);

r1 = r1 + 100;
th1 = 90 - (180./3.1415926*th1);
ph1 = 180./3.1415926*ph1;



% Make 3-d density arrays
index = 1
nthp = nth+1
for k = 1:nr
   for j = 1:nth
      for i = 1:nph
	  ked(i,j,k) = kedep(index);
          index = index + 1;
      end
   end
end

      ked = log10(ked);

for k = 1:nr
   for j = 1:nth
      ked(nph,j,k) = ked(nph-1,j,k);
      ked(1,j,k) = ked(2,j,k);
%ked(nph,j,k) = ked(2,j,k);
   end
end





% Use the slice routine to make the contours

lonlatplane = 140.0

	g = slice(th1,ph1,r1,ked, 0, 0, lonlatplane );
%	g = slice(th1,ph1,r1,ne, [], [], [200 400 600 800 1000 1200 1400] );
%	g = slice(ph1,th1,r1,ne, [], [], lonlatplane );
%g = slice(th1,ph1,r1,ked,90-(180./3.1415926*(thb1:dth:(nth)*dth)) ,...
%180./3.1415926*(0:dph:(nph-1)*dph), 100);
%g = slice(th1,ph1,r1,ked,90-(180./3.1415926*(thb1:dth:(nth)*dth)) ,...
%0, 100);
%g = slice(th1,ph1,r1,ked,0,0, 100+(0:dr:(nr-1)*dr));

set(0,'defaultaxesfontsize',20);
set(g,'FaceColor','interp','facelighting','none','EdgeColor','none',...
      'DiffuseStrength',0.0,'AmbientStrength',1,...
    'SpecularColorReflectance',0,'SpecularStrength',0.0)
set(gca,'XDir','rev')

xlabel('Latitude (deg.)')
ylabel('Longitude (deg.)')
zlabel('Altitude (km)')
%title({'Ion Energy Deposition' 'Log(eV cm^-^3 s^-^1), Alt. = 110 km'})
%title({'SW H^+ Energy Deposition' 'Log(eV cm^-^3 s^-^1), Alt. = 195 km'})
title({'Ion Energy Deposition' 'Log(eV cm^-^3 s^-^1), Alt. = 100 km'})


view(-90,90)                     % Long. vs. Latitude
view(0,0)                        % Lat. vs. Altitude
axis([-90 90 0 360 100 450])


colorbar
caxis([0 5.2])
daspect([1 1 1])

% Color of the axis labels

%set(gca,'Color',[0.1 0.1 0.1]);      % plot background to black

%set(gca,'XColor',[0.847059 0.160784 0])    % Plot labels to red
%set(gca,'YColor',[0.847059 0.160784 0])    % Plot labels to red
%set(gca,'ZColor',[0.847059 0.160784 0])    % Plot labels to red

view(90,90)

caxis([0 4.5])


stop

figure

ii = fix((lonlatplane - 100)/dr)

nosph = zeros(nph,nth);
no2sph = zeros(nph,nth);
nesph = zeros(nph,nth);

nthp = nth + 1;

for j = 1:nth;
   lat0 = (j-1)*dth;
   lat0 = 90 - (180./3.1415926*lat0);
   for k = 1:nph
	  lon0 = (k-1)*dph * (180./3.1415926);
          lat(k,j) = lat0;
          lon(k,j) = lon0;
          nosph(k,j)  = nop(k,j,ii);
          no2sph(k,j) = no2p(k,j,ii);
	  nesph(k,j)  = ne(k,j,ii);
   end
end
%for j = 1:nth
%   for k = 1:nph
%          nosph(k,j)  = nop(k,nthp-j,ii);
%          no2sph(k,j) = no2p(k,nthp-j,ii);
%	  nesph(k,j)  = ne(k,nthp-j,ii);
%   end
%end


%lat0 = 0:dth:(nth-1)*dth;
 %     lon0 = 0:dph:(nph-1)*dph;

%[lat,lon] = meshgrid(lat0,lon0,1);
%lat = 90 - (180./3.1415926*th1);
%lon = 180./3.1415926*ph1;


surf(lon,lat,nesph);

xlabel('Latitude (deg.)')
ylabel('Longitude (deg.)')
zlabel('Log [ Ionospheric Density (cm^-^3) ]')
title({'Ionospheric Density' 'Altitude = 1000 km'})

axis([0 360 -89 89 -5 5])
s.EdgeColor = 'none';    
colorbar
shading interp
lighting phong