% Dibuixo el potencial electrostatic a partir d'un arxiu netCDF procedent d'AMIE.
clear all
close all
k = 0;
ncid = netcdf.open('/hao/aim2/marsal/GangLu/amie/sep16_2005_nh.nc', '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
FAC = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'FAC')); % Llegeixo la matriu pot i li assigno la variable FAC. pot 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.
lon0 = pi/180*netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lon'))';
lon = [lon0, lon0(1) + 2*pi]; % Afegeixo la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
colat = pi/180*(90 - netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lat'))); % Aqui no cal ivertir el vector latituds, ja que comenca al Pol Nord.
[lonMat, colatMat] = meshgrid(lon, colat); % Genero la malla.
PosUT0 = 1: 20: 288;
for PosUT = PosUT0
    %PosUT = 217; % Posicio dintre del vector ut.
    %Posilev = 1; % A AMIE no hi ha nivells d'alcada.
    UT = UTvec(PosUT); % UT es ara el temps universal pel qual es fara el dibuix.
    %FAC = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'pot')); % Llegeixo la matriu pot i li assigno la variable FAC. pot 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.
    %colat = colat(1: end/2); El rang de latituds amb AMIE es mes reduit (de 50 lat al Pol).
    FAC_t0 = NaN(size(FAC, 1), size(FAC, 2));
    FAC_t0(:) = FAC(:, :, PosUT); % M'interessa el potencial per a un temps UT determinat.
    FAC_t0 = FAC_t0'; % Aqui no em cal invertir amb flipud, ja que la latitud comenca al Pol Nord.
    %FAC_t = flipud(FAC_t'); % Inverteixo amb flipud, ja que la latitud comen?a a l'hemisferi Sud.
    %FAC_t = FAC_t(1: end/2, :); % M'interessa nom?s l'hemisferi Nord.
    FAC_t = [FAC_t0, FAC_t0(:, 1)]; % Afegeixo la FAC corresponent a la longitud positiva coincident amb la longitud inicial + 360 graus. De cara a dibuixar.

    % Dibuixo les dades de FAC en un mapa polar:
    %figure
    [x_0, y_0] = pol2cart(lonMat + 15*pi/180*UT, 180/pi*colatMat); % Passo les matrius theta i lambda a cartesianes. x_0 i y_0 son matrius 91x25 (de 0 a 90 graus i de -180 a 180 graus).
    h = polar(-pi: pi/100: pi, 40*ones(1, 201), '-k'); delete(h)
    hold on
    contourf(x_0, y_0, FAC_t, 'LineStyle', 'none', 'DisplayName', 'FAC(V)')

    % Pinto meridians i paral.lels geografics:
    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
    %polar((24-UT)/12*pi*[1 1], [0 40], 'k') % Marco amb una linia fosca el meridia de mitjanit.
    xlim([-40 40]); ylim([-40 40])
    xlabel('18 LT', 'FontSize', 8);
    ylabel('0 LT', 'FontSize', 8);
    %view(90 - 15*UT, 90) % En principi era per girar el mapa amb la finalitat que mitjanit quedi a baix, pero no funciona be!
    view(90, 90)
    axis equal
    k = k + 1;
    M(k) = getframe;
end
movie(M, 1, 3)