% SIMPLECLOUD_PLOT   Automatic plotting of cloud scenario.
%
%    A function to visualize a cloud case generated by SIMPLECLOUD_INIT.
%    The function adapts automatically to the dimesnion, position and
%    size of the cloud.
%
%    Note that the cloud is plotted in latitude and longitude, and not
%    in geoemtrical distance. This means that, for example, a cloud with 
%    quadratic horisontal extension can look rectangular in the view from
%    above.
%     
% FORMAT   simplecloud_plot(C)
%
% IN       C      Settings from SIMPLECLOUD_INIT.

% HISTORY: 2003-03-09  Created by Patrick Eriksson


function simplecloud_plot(C)


%=== Get an upper limit for plotting
%=== Plots are from 0 km to highest altitude with cloud/rain*WFACTOR
%
wfactor = 1.2;   
%
zmax = -1;
%
if C.do_iwc 
  zmax = max( C.iwc(:,1) );
end
%
if C.do_lwc 
  zmax = max( [ zmax max( C.lwc(:,1) ) ] );
end
%
if C.do_rain
  zmax = max( [ zmax max( C.rain(:,1) ) ] );
end



%=== Create some help variables
%
z = 0 : 100: zmax*1.2;
%
if C.dim == 1
  c = zeros( length(z), 3 );
  for i = 1 : length(z)
    [iwc,lwc,rain] = simplecloud_get( C, z(i) );
    c(i,:) = [ iwc, lwc, rain ];
  end
else
  if C.hshape == 1
    lat = C.lat + linspace( wfactor*(-C.latlength/2-C.hedge), ...
                                wfactor*(C.latlength/2+C.hedge), 100 ) / 111e3;
    lon = C.lon + linspace( wfactor*(-C.lonlength/2-C.hedge), ...
                wfactor*(C.lonlength/2+C.hedge), 100 )/111e3/cos(C.lat*pi/180);
  else
    lat = C.lat + linspace( -C.latlength*1.5, C.latlength*1.5, 100 ) / 111e3;
    lon = C.lon + linspace( -C.lonlength*1.5, C.lonlength*1.5, 100 ) / ...
                                                       111e3/cos(C.lat*pi/180);
  end

  %- Latitude cross-section
  clat = zeros( length(z), length(lat), 3 );
  for i = 1 : length(z)
    for j = 1 : length(lat)
      [iwc,lwc,rain] = simplecloud_get( C, z(i), lat(j), C.lon );
      clat(i,j,:) = [ iwc, lwc, rain ];
    end
  end
  %- Longitude cross-section
  clon = zeros( length(z), length(lon), 3 );
  for i = 1 : length(z)
    for j = 1 : length(lon)
      [iwc,lwc,rain] = simplecloud_get( C, z(i), C.lat, lon(j) );
      clon(i,j,:) = [ iwc, lwc, rain ];
    end
  end
end


fig = gcf - 1;



%=== IWC
%
if C.do_iwc
    if C.dim == 1        
        %
        fig = fig + 1; figure( fig ); 
        plot( c(:,1), z/1e3 )
        xlabel('Ice water content [g/m3]', 'FontSize', 12 )
        ylabel('Altitude [km]', 'FontSize', 12 )
        axis([0 max(C.iwc(:,2))*1.1 0 zmax*wfactor/1e3])
        %        
    else
        %
        fig = fig + 1; figure( fig ); 
        contourf( lat, z/1e3, squeeze( clat(:,:,1) ) );
        shading flat
        colorbar
        axis([lat(1) tail(lat,1) 0 zmax*wfactor/1e3])
        xlabel('Latitude [deg]', 'FontSize', 12);
        ylabel('Altitude [km]', 'FontSize', 12 );
        title( 'Latitude cross-section of ice water content [g/m3]', ...
            'FontSize', 12 )
        %
        fig = fig + 1; figure( fig ); 
        contourf( lon, z/1e3, squeeze( clon(:,:,1) ) );
        shading flat
        colorbar
        axis([lon(1) tail(lon,1) 0 zmax*wfactor/1e3])
        xlabel('Longitude [deg]', 'FontSize', 12);
        ylabel('Altitude [km]', 'FontSize', 12 );
        title( 'Longitude cross-section of ice water content [g/m3]', ...
            'FontSize', 12 )
        %
        [u,iz] = max( C.iwc(:,2) );
        cz = zeros( length(lat), length(lon) );
        for i = 1 : length(lat)
            for j = 1 : length(lon)
                [iwc,lwc,rain] = simplecloud_get( C, C.iwc(iz,1), lat(i), lon(j) );
                cz(i,j) = iwc;
            end
        end
        fig = fig + 1; figure( fig ); 
        contourf( lon, lat, cz );
        shading flat
        colorbar
        axis([ lon(1) tail(lon,1) lat(1) tail(lat,1) ])
        axis equal
        ylabel('Latitude [deg]', 'FontSize', 12);
        xlabel('Longitude [deg]', 'FontSize', 12);
        title( sprintf('Ice water content at %.2f km [g/m3]', C.iwc(iz,1)/1e3 ), ...
            'FontSize', 12 )
        
    end
end



%=== LWC
%
if C.do_lwc
    if C.dim == 1
        %
        fig = fig + 1; figure( fig ); 
        plot( c(:,2), z/1e3 )
        xlabel('Liquid water content [g/m3]', 'FontSize', 12 )
        ylabel('Altitude [km]', 'FontSize', 12 )
        axis([0 max(C.lwc(:,2))*1.1 0 zmax*wfactor/1e3])
        %
    else
        %
        fig = fig + 1; figure( fig ); 
        contourf( lat, z/1e3, squeeze( clat(:,:,2) ) );
        shading flat
        colorbar
        axis([lat(1) tail(lat,1) 0 zmax*wfactor/1e3])
        xlabel('Latitude [deg]', 'FontSize', 12);
        ylabel('Altitude [km]', 'FontSize', 12 );
        title( 'Latitude cross-section of liquid water content [g/m3]', ...
            'FontSize', 12 )
        %
        fig = fig + 1; figure( fig ); 
        contourf( lon, z/1e3, squeeze( clon(:,:,2) ) );
        shading flat
        colorbar
        axis([lon(1) tail(lon,1) 0 zmax*wfactor/1e3])
        xlabel('Longitude [deg]', 'FontSize', 12);
        ylabel('Altitude [km]', 'FontSize', 12 );
        title( 'Longitude cross-section of liquid water content [g/m3]', ...
            'FontSize', 12 )
        %
        [u,iz] = max( C.lwc(:,2) );
        cz = zeros( length(lat), length(lon) );
        for i = 1 : length(lat)
            for j = 1 : length(lon)
                [iwc,lwc,rain] = simplecloud_get( C, C.lwc(iz,1), lat(i), lon(j) );
                cz(i,j) = lwc;
            end
        end
        fig = fig + 1; figure( fig ); 
        contourf( lon, lat, cz );
        shading flat
        colorbar
        axis([ lon(1) tail(lon,1) lat(1) tail(lat,1) ])
        axis equal
        ylabel('Latitude [deg]', 'FontSize', 12);
        xlabel('Longitude [deg]', 'FontSize', 12);
        title( ...
            sprintf('Liquid water content at %.2f km [g/m3]', C.lwc(iz,1)/1e3 ), ...
            'FontSize', 12 )
    end
end



%=== RAIN
%
if C.do_rain
    if C.dim == 1
        %
        fig = fig + 1; figure( fig ); 
        plot( c(:,3), z/1e3 )
        xlabel('Rain rate [mm/h]', 'FontSize', 12 )
        ylabel('Altitude [km]', 'FontSize', 12 )
        axis([0 max(C.rain(:,2))*1.1 0 zmax*wfactor/1e3])
        %
    else
        %
        fig = fig + 1; figure( fig ); 
        contourf( lat, z/1e3, squeeze( clat(:,:,3) ) );
        shading flat
        colorbar
        axis([lat(1) tail(lat,1) 0 zmax*wfactor/1e3])
        xlabel('Latitude [deg]', 'FontSize', 12);
        ylabel('Altitude [km]', 'FontSize', 12 );
        title( 'Latitude cross-section of rain rate [mm/h]', ...
            'FontSize', 12 )
        %
        fig = fig + 1; figure( fig ); 
        contourf( lon, z/1e3, squeeze( clon(:,:,3) ) );
        shading flat
        colorbar
        axis([lon(1) tail(lon,1) 0 zmax*wfactor/1e3])
        xlabel('Longitude [deg]', 'FontSize', 12);
        ylabel('Altitude [km]', 'FontSize', 12 );
        title( 'Longitude cross-section of rain rate [mm/h]', ...
            'FontSize', 12 )
        %
        [u,iz] = max( C.rain(:,2) );
        cz = zeros( length(lat), length(lon) );
        for i = 1 : length(lat)
            for j = 1 : length(lon)
                [iwc,lwc,rain] = simplecloud_get( C, C.rain(iz,1), lat(i), lon(j) );
                cz(i,j) = rain;
            end
        end
        fig = fig + 1; figure( fig ); 
        contourf( lon, lat, cz );
        shading flat
        colorbar
        axis([ lon(1) tail(lon,1) lat(1) tail(lat,1) ])
        axis equal
        ylabel('Latitude [deg]', 'FontSize', 12);
        xlabel('Longitude [deg]', 'FontSize', 12);
        title( ...
            sprintf('Rain rate at %.2f km [mm/h]', C.rain(iz,1)/1e3 ), ...
            'FontSize', 12 )
    end
end