% Dibuixo el potencial electrostatic a partir d'un arxiu netCDF procedent dels TGCMs. Faig una pel.licula per a diferents temps.

clear all
close all

% Parametres d'entrada a modificar per l'usuari:
Posilev = 29; % Posicio dintre del vector ilev (dona l'alcada com un nivell de pressio constant dintre de la ionosfera).
PosUTvec = 1: 210; % Posicions dintre del vector UT per a les que es dibuixara el potencial.
Field = 'HALL'; % Camp a dibuixar.
Unitats = 'mho'; % Unitats del camp a dibuixar.
LimitsFLD = [-5, 30]; % Limits de l'escala de colors del dibuix.
nomsarxnc{1} = './tiegcm/s_sep05_fism_093.nc'; % Nom de l'arxiu netCDF (*.nc) a obrir.
nomsarxnc{2} = './tiegcm/s_sep05_fism_094.nc'; % Nom de l'arxiu netCDF (*.nc) a obrir.
nomsarxnc{3} = './tiegcm/s_sep05_fism_095.nc'; % Nom de l'arxiu netCDF (*.nc) a obrir.
nomsarxnc{4} = './tiegcm/s_sep05_fism_096.nc'; % Nom de l'arxiu netCDF (*.nc) a obrir.
nomsarxnc{5} = './tiegcm/s_sep05_fism_097.nc'; % Nom de l'arxiu netCDF (*.nc) a obrir.
nomsarxnc{6} = './tiegcm/s_sep05_fism_098.nc'; % Nom de l'arxiu netCDF (*.nc) a obrir.
nomsarxnc{7} = './tiegcm/s_sep05_fism_099.nc'; % Nom de l'arxiu netCDF (*.nc) a obrir.
dimlon = 'lon'; % dimensio longitud utilitzada pel camp Field. Cal fer un ncdump -c nomarxnc per veure quina dimensio utilitza!
dimlat = 'lat'; % dimensio latitud utilitzada pel camp Field. Cal fer un ncdump -c nomarxnc per veure quina dimensio utilitza!
MinColatProj = 40; % Minima colatitud de la projeccio, en graus.

% Obro l'arxiu netCDF, llegeixo el que m'interessa i tanco:
ncid = netcdf.open(nomsarxnc{1}, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc)
IdVarField = netcdf.inqVarID(ncid, Field);
[varname, xtype, dimids, natts] = netcdf.inqVar(ncid, IdVarField);
disp(['Nom i longitud de cada dimensio de la variable ', Field, ':']) 
for dimid = dimids
    [dimname, dimlen] = netcdf.inqDim(ncid, dimid);
    disp([dimname, ', ', num2str(dimlen)])
end
UTvec = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ut')); % Llegeixo el vector ut i li assigno la variable UT
FLD = netcdf.getVar(ncid, netcdf.inqVarID(ncid, Field)); % Llegeixo la variable especificada i li assigno la variable FLD.
lon_coor = pi/180*netcdf.getVar(ncid, netcdf.inqVarID(ncid, dimlon))'; % longitud utilitzada en el camp Field (pot ser geografica, magnetica, etc.).
%lon_coor = pi/180*netcdf.getVar(ncid, 1)'; % Podria automatitzar-ho aixi, pero m'estimo mes assegurar-me fent un ncdump i veient de quines dimensions depen la variable Field.
lon_coor = [lon_coor, lon_coor(1) + 2*pi]; % Afegeixo la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
colat_coor = pi/180*(90 - wrev(netcdf.getVar(ncid, netcdf.inqVarID(ncid, dimlat)))); % latitud utilitzada en el camp Field (pot ser geografica, magnetica, etc.). Cal invertir el vector de latituds pr?viament, ja que comen?a a l'hemisferi Sud.
netcdf.close(ncid)

% Llegeixo la resta d'arxius i redefineixo les matrius UTvec i FLD.
for numarxnc = 2: length(nomsarxnc)
    ncid = netcdf.open(nomsarxnc{numarxnc}, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc)
    UTvec = cat(1, UTvec, netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ut'))); % Llegeixo el vector ut i li assigno la variable UT
    FLD = cat(4, FLD, netcdf.getVar(ncid, netcdf.inqVarID(ncid, Field))); % Llegeixo la variable especificada i li assigno la variable FLD.
    netcdf.close(ncid)
end

% M'arreglo el que he llegit a l'arxiu segons els meus interessos:
UT = UTvec(PosUTvec(1)); % UT ?s ara el temps universal pel qual es far? el dibuix.
colat_coor = colat_coor(1: end/2); % M'interessa nomes l'hemisferi Nord.
colat_coor(colat_coor > pi/180*MinColatProj) = NaN;
[lonMat, colatMat] = meshgrid(lon_coor, colat_coor); % Genero la malla.
FLD_t_ilev0 = NaN(size(FLD, 1), size(FLD, 2));
FLD_t_ilev0(:) = FLD(:, :, Posilev, PosUTvec(1)); % M'interessa el potencial per a una al?ada determinada (donada per ilev) i per a un temps UT determinat.
FLD_t_ilev1 = flipud(FLD_t_ilev0'); % Inverteixo amb flipud, ja que la latitud comen?a a l'hemisferi Sud.
FLD_t_ilev2 = FLD_t_ilev1(1: end/2, :); % M'interessa nom?s l'hemisferi Nord.
FLD_t_ilev = [FLD_t_ilev2, FLD_t_ilev2(:, 1)]; % Afegeixo FLD corresponent a la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
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([-MinColatProj MinColatProj])
ylim([-MinColatProj MinColatProj])
xlabel('0 LT', 'FontSize', 8);
ylabel('18 LT', 'FontSize', 8);
colorbar
hold on

% Passo les matrius theta i lambda a cartesianes:
[x_0, y_0] = pol2cart(lonMat + 15*pi/180*UT - pi/2, 180/pi*colatMat); % 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_ilev, 'LineStyle', 'none', 'DisplayName', [Field, '(', Unitats, ')'], 'LevelStep', 1);
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_ilev')

% Afegeixo linies de costa:
load coast
lonmap = pi/180*long; % rad
colatmap = pi/180*(90 - lat); % rad
clear long lat
colatmap(colatmap > pi/180*MinColatProj) = NaN;
[x_m, y_m] = pol2cart(lonmap + 15*pi/180*UT - pi/2, 180/pi*colatmap); % 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?).
hmap = plot(x_m, y_m, 'k');
set(hmap, 'XDataSource', 'x_m', 'YDataSource', 'y_m')

% 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 MinColatProj],  ':k')
end

% Itero per a cada pas de temps:
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_t_ilev0(:) = FLD(:, :, Posilev, PosUT); % M'interessa el potencial per a una al?ada determinada (donada per ilev) i per a un temps UT determinat.
    FLD_t_ilev1 = flipud(FLD_t_ilev0'); % Inverteixo amb flipud, ja que la latitud comen?a a l'hemisferi Sud.
    FLD_t_ilev2 = FLD_t_ilev1(1: end/2, :); % M'interessa nom?s l'hemisferi Nord.
    FLD_t_ilev = [FLD_t_ilev2, FLD_t_ilev2(:, 1)]; % Afegeixo la FLD corresponent a la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
    
    [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_m, y_m] = pol2cart(lonmap + 15*pi/180*UT - pi/2, 180/pi*colatmap); % 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

% Presento la pel.licula:
movie(MT, 1, 5)