% Dibuixo el potencial electrost?tic a partir d'un arxiu netCDF procedent dels TGCMs. Faig una pel?l?cula per a diferents temps.
clear all
close all
% Par?metres:
%Posilev = 29; % Posici? dintre del vector ilev (d?na l'al?ada com un nivell de pressi? constant dintre de la ionosfera).
PosUTvec = 19: 19; % Posici? dintre del vector UT per al que es dibuixar? el potencial.
Field = 'POTEN';
Unitats = 'V';
LimitsFLD = [-5e4, 5e4];
%nomarxnc = '/hao/aim2/marsal/tiegcm/tiegcm-linux/MARSAL.tiegcm_trunk.ssep1605.nc';
nomarxnc = '/hao/aim2/marsal/tiegcm/tiegcm-linux/TGCM.tiegcm1.93.pclim_smin_011.nc';

ncid = netcdf.open(nomarxnc, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc)
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 variable especificada i li assigno la variable FLD.
mlon = pi/180*netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'mlon'))';
%mlon = [mlon, mlon(1) + 2*pi]; % Afegeixo la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
mcolat = pi/180*(90 - wrev(netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'mlat')))); % Cal invertir el vector de latituds pr?viament, ja que comen?a a l'hemisferi Sud.
mcolat = mcolat(1: ceil(end/2)); % M'interessa nom?s l'hemisferi Nord.
[mlonMat, mcolatMat] = meshgrid(mlon, mcolat); % Genero la malla.
FLD_t0 = NaN(size(FLD, 1), size(FLD, 2));
FLD_t0(:) = FLD(:, :, PosUTvec(1)); % M'interessa el potencial per a una al?ada determinada (donada per ilev) i per a un temps UT determinat.
FLD_t1 = flipud(FLD_t0'); % Inverteixo amb flipud, ja que la latitud comen?a a l'hemisferi Sud.
FLD_t2 = FLD_t1(1: ceil(end/2), :); % M'interessa nom?s l'hemisferi Nord.
%FLD_t = [FLD_t2, FLD_t2(:, 1)]; % Afegeixo FLD corresponent a la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
FLD_t = FLD_t2;
MT(1, length(PosUTvec)) = struct('cdata', [], 'colormap', []); % Inicialitzo l'estructura MT.

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

[x_0, y_0] = pol2cart(mlonMat + 15*pi/180*UT - pi/2, 180/pi*mcolatMat); % 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 FLD en un mapa polar:
[C, hc] = contourf(x_0, y_0, FLD_t, 'LineStyle', 'none', 'DisplayName', [Field, '(', Unitats, ')'], 'LevelStep', .02);
set(gca, 'CLim', LimitsFLD) % 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, 20*ones(1, 201), ':k'); polar(-pi: pi/100: pi, 40*ones(1, 201), ':k');
for i = 0: 30: 359
    polar([i i]*pi/180, [0 40],  ':k')
end
MT(1) = getframe;
k = 1;
for PosUT = PosUTvec(2: end)
    k = k + 1;
    UT = UTvec(PosUT); % UT ?s ara el temps universal pel qual es far? el dibuix.
    FLD_t0(:) = FLD(:, :, PosUT); % M'interessa el potencial per a una al?ada determinada (donada per ilev) i per a un temps UT determinat.
    FLD_t1 = flipud(FLD_t0'); % Inverteixo amb flipud, ja que la latitud comen?a a l'hemisferi Sud.
    FLD_t2 = FLD_t1(1: ceil(end/2), :); % M'interessa nom?s l'hemisferi Nord.
    %FLD_t = [FLD_t2, FLD_t2(:, 1)]; % Afegeixo la FLD corresponent a la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
    FLD_t = FLD_t2;
    
    [x_0, y_0] = pol2cart(mlonMat + 15*pi/180*UT - pi/2, 180/pi*mcolatMat); % 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
    MT(k) = getframe;
end
movie(MT, 1, 5)