% Dibuixo el FLDencial electrostatic a partir d'un arxiu netCDF procedent d'AMIE.
clear all
close all
% Par?metres:
PosUTvec = 1: 8; % Posicio dintre del vector ut.
Field = 'magpot_tot';
Unitats = 'nT.m';
%LlistaNivells = -0.12: 0.02: 0.12;
nomarxnc = '/home/marsal/maggrd/maggrd.d/mag_timegcm_doy339_solF_qnm.nc';

ncid = netcdf.open(nomarxnc, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc). nh -> hemisferi nord.
%UTvec = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ut')); % Llegeixo el vector ut i li assigno la variable UT
%UT = UTvec(PosUTvec(1)); % UT ?s ara el temps universal pel qual es far? el dibuix.
FLD = netcdf.getVar(ncid, netcdf.inqVarID(ncid, Field)); % Llegeixo la matriu FLD i li assigno la variable FLD. FLD es aqui una matriu 3D (lonxlatxtime). Compte! sembla que esta a l'inreves a l'especificacio de l'arxiu netCDF quan es fa un "ncdump -c sep16_2005_nh.nc" al linux.
rmlt = pi/12*netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'rmlt'))'; % En realitat ?s MLT en caomptes de longitud!!!
%lon = [lon, lon(1) + 2*pi]; % Afegeixo la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
rmcolat = pi/180*(90 - wrev(netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'rmlat')))); % Aqui no cal ivertir el vector latituds, ja que comenca al Pol Nord.
rmcolat = rmcolat(1: ceil(length(rmcolat)/2)); % M'interessa nom?s l'hemisferi Nord.
[rmltMat, rmcolatMat] = meshgrid(rmlt, rmcolat); % Genero la malla.
FLD_t0 = NaN(size(FLD, 1), size(FLD, 2));
FLD_t0(:) = FLD(:, :, PosUTvec(1)); % M'interessa el FLDencial per a un temps UT determinat.
FLD_t1 = flipud(FLD_t0'); % Aqui no em cal invertir amb flipud, ja que la latitud comenca al Pol Nord.
FLD_t2 = FLD_t1(1: ceil(size(FLD_t1, 1)/2), :); % Afegeixo la Phi corresponent a la longitud positiva coincident amb la longitud inicial + 360 graus. De cara a dibuixar.
FLD_t = FLD_t2;
MM(1, length(PosUTvec)) = struct('cdata', [], 'colormap', []); % Inicialitzo l'estructura MA.

% Estableixo els nivells de colors adequats i els seus limits:
minFLD = min(FLD(:)); maxFLD = max(FLD(:));
Salt = (maxFLD - minFLD)/12; Salt = round(Salt*10^(-floor(log10(Salt))))*10^(floor(log10(Salt)));
LlistaNivellsMin = Salt*floor(minFLD/Salt); LlistaNivellsMax = Salt*ceil(maxFLD/Salt);
LlistaNivells = LlistaNivellsMin: Salt: LlistaNivellsMax;

% Creo un mapa polar:
h = polar(-pi: pi/100: pi, 30*ones(1, 201), '-k'); delete(h)
xlim([-90 90]);
ylim([-90 90])
xlabel('0 MLT', 'FontSize', 8);
ylabel('18 MLT', 'FontSize', 8);
colorbar
hold on

%[x_0, y_0] = pol2cart(lonMat + 15*pi/180*UT - pi/2, 180/pi*colatMat); % Passo les matrius theta i lambda a cartesianes. x_0 i y_0 s?n matrius 91x25 (de 0? a 90? i de -180? a 180?).
[x_0, y_0] = pol2cart(rmltMat - pi/2, 180/pi*rmcolatMat); % Passo les matrius theta i lambda a cartesianes. x_0 i y_0 s?n matrius 91x25 (de 0? a 90? i de -180? a 180?).
% Dibuixo les dades de Phi en un mapa polar:
[C, hc] = contourf(x_0, y_0, FLD_t, 'LineStyle', 'none', 'DisplayName', [Field, '(', Unitats, ')'], 'LevelList', LlistaNivells);
set(gca, 'CLim', [LlistaNivells(1), LlistaNivells(end)]) % Estableixo els l?mits de colors de l'escala vertical a fi que siguin iguals per als diferents passos de temps.
set(hc, 'XDataSource', 'x_0', 'YDataSource', 'y_0', 'ZDataSource', 'FLD_t')

% Pinto meridians i paral?lels geogr?fics:
polar(-pi: pi/100: pi, 30*ones(1, 201), ':k'); polar(-pi: pi/100: pi, 60*ones(1, 201), ':k');
for i = 0: 30: 359
    polar([i i]*pi/180, [0 90],  ':k')
end
MM(1) = getframe;
k = 1;
for PosUT = PosUTvec(2: end)
    %UT = UTvec(PosUT); % UT es ara el temps universal pel qual es fara el dibuix.
    FLD_t0(:) = FLD(:, :, PosUT); % M'interessa el FLDencial per a un temps UT determinat.
    FLD_t1 = flipud(FLD_t0'); % Aqui no em cal invertir amb flipud, ja que la latitud comenca al Pol Nord.
    FLD_t2 = FLD_t1(1: ceil(size(FLD_t1, 1)/2), :); % Afegeixo la Phi corresponent a la longitud positiva coincident amb la longitud inicial + 360 graus. De cara a dibuixar.
    FLD_t = FLD_t2;

    % Dibuixo les dades de Phi en un mapa polar:
%    figure
%    [x_0, y_0] = pol2cart(lonMat + 15*pi/180*UT - pi/2, 180/pi*colatMat); % Passo les matrius theta i lambda a cartesianes. x_0 i y_0 s?n matrius 91x25 (de 0? a 90? i de -180? a 180?).
    [x_0, y_0] = pol2cart(rmltMat - pi/2, 180/pi*rmcolatMat); % Passo les matrius theta i lambda a cartesianes. x_0 i y_0 s?n matrius 91x25 (de 0? a 90? i de -180? a 180?).
    refreshdata
%    drawnow
    k = k + 1;
    MM(k) = getframe;
end
movie(MM, 1, 5)