%
% Dibuixo el camp Field a partir d'un arxiu netCDF procedent dels TGCMs. Faig una pel·lícula per a diferents temps.
%
clear all
close all

% Paràmetres d'entrada a modificar per l'usuari:
PosUTvec = 1: 15: 720; % Posició dintre del vector UT per al que es dibuixarà el potencial.
Field = 'jr';
Field_corr = 'FAC';
Unitats = 'A/m2';
any = 2010.7;
DOY = 237;
nomarxnc = '.\20100825.north.nc';
%dimlon = 'lon'; % Dimensió longitud utilitzada pel camp Field. Cal fer un 'ncdump -c nomarxnc' (a Linux) per veure quina dimensió utilitza la variable Field!
%dimlat = 'lat'; % Dimensió latitud utilitzada pel camp Field. Cal fer un 'ncdump -c nomarxnc' (a Linux) per veure quina dimensió utilitza la variable Field!
MaxColatProj = 60; % º. Mínima colatitud de la projecció.

% Obro l'arxiu netCDF, llegeixo les variables que m'interessen i tanco:
ncid = netcdf.open(nomarxnc, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc)
IdVarField = netcdf.inqVarID(ncid, Field); % Retorna la IDentitat de la variable Field.
[varname, xtype, dimids, natts] = netcdf.inqVar(ncid, IdVarField); % M'interessen les identitats de les diferents dimensions de la variable Field.
disp(['Nom i longitud de cada dimensio (o argument) de la variable ', Field, ':'])
for dimid = dimids
    [dimnom, dimllarg] = netcdf.inqDim(ncid, dimid);
    disp([dimnom, ', ', num2str(dimllarg)])
end
UTvec = single(netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'start_hr'))) + single(netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'start_mt')))/60 + single(netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'start_sc')))/3600; % Llegeixo el vector ut i li assigno la variable UTvec.
FLD = 1e-6*netcdf.getVar(ncid, netcdf.inqVarID(ncid, Field)); % Llegeixo la variable especificada i li assigno la variable FLD.
nlon = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'nlon')); nlon = nlon(1);
nlat = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'nlat')); nlat = nlat(1);
if nlon*nlat*length(UTvec) ~= numel(FLD)
    Alguna cosa no quadra
end
mlt0 = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'mlt')); % Longitud utilitzada en el camp Field (pot ser geografica, magnetica, mlt, etc.).
mlt0(mlt0 == 24) = 0;
colat0 = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'colat')); % Longitud utilitzada en el camp Field (pot ser geografica, magnetica, mlt, etc.).
colatN = unique(colat0);
sinI_N = repmat(2*cosd(colatN)./sqrt(1 + 3*cosd(colatN).^2), 1, nlon);
netcdf.close(ncid)

% 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.
lonMat = NaN(nlat, nlon);
lonMat(:) = pi/12*mlt0(:, 1);
%lonMat = pi/180*netcdf.getVar(ncid, 1)'; % Podria automatitzar-ho així, però m'estimo més assegurar-me fent un ncdump i veient de quines dimensions depen la variable Field.
lonMat = [lonMat, lonMat(:, 1) + 2*pi]; % Afegeixo la longitud positiva coincident amb la longitud inicial + 360º. De cara a dibuixar.
colatMat = NaN(nlat, nlon);
colatMat(:) = pi/180*colat0(:, 1);
colatMat = [colatMat, colatMat(:, 1)];
FLD_t0 = NaN(nlat, nlon);
FLD_t0(:) = FLD(:, PosUTvec(1));
FLD_t0 = FLD_t0./abs(sinI_N + eps); % Desfaig la projecció radial. Ara tinc FACs.
FLD_t = double([FLD_t0, FLD_t0(:, 1)]); % Afegeixo FLD corresponent a la longitud positiva coincident amb la longitud inicial + 360º. De cara a dibuixar.

MAM(1, length(PosUTvec)) = struct('cdata', [], 'colormap', []); % Inicialitzo l'estructura MAM.

% 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;
%LimitsFLD = [LlistaNivellsMin, LlistaNivellsMax];
LimitsFLD = [-1.4e-6, .8e-6];
Salt = .2e-6;
LlistaNivells = LimitsFLD(1): Salt: LimitsFLD(2);

% Creo un mapa polar:
h = polar(-pi: pi/100: pi, 30*ones(1, 201), '-k'); delete(h)
xlim([-MaxColatProj MaxColatProj])
ylim([-MaxColatProj MaxColatProj])
xlabel('00 MLT', 'FontSize', 8);
ylabel('18 MLT', '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º).
colatMat(colatMat > pi/180*MaxColatProj) = NaN;
[x_0, y_0] = pol2cart(lonMat - 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, 'LineStyle', 'none', 'DisplayName', [Field_corr, ' (', Unitats, ')'], 'LevelList', LlistaNivells);
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')

% Afegeixo línies de costa Compte! El mapa està en un sistema de coordenades geomagnètiques; no Apex.
load coast
[lonGmMap latGmMap lonPolGm colatPolGm] = Geog2Geom(any, long, lat); % º. Converteixo la malla geogràfica a geomagnètica.
clear long lat
lonGmMap = lonGmMap*pi/180; %rad
colatGmMap = pi/180*(90 - latGmMap); %rad
colatGmMap(colatGmMap > pi/180*MaxColatProj) = NaN;
%MUT = UT + lonPolGm/15;
MUT = MLT(lonPolGm, 0, any, DOY, UT, 1);
[x_m, y_m] = pol2cart(lonGmMap + pi/12*MUT - pi/2, 180/pi*colatGmMap);
hmap = plot(x_m, y_m, 'b');
set(hmap, 'XDataSource', 'x_m', 'YDataSource', 'y_m')

% Pinto meridians i paral·lels geogràfics:
for i = 0: 30: 359
    polar([i i]*pi/180, [0 MaxColatProj],  ':k')
end
polar(-pi: pi/100: pi, 20*ones(1, 201), ':k'); polar(-pi: pi/100: pi, 40*ones(1, 201), ':k');

% Itero per a cada pas de temps:
MAM(1) = getframe; % Agafa la imatge i la guarda en memòria.
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.
    %MUT = UT + lonPolGm/15;
    MUT = MLT(lonPolGm, 0, any, DOY, UT, 1);
    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_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 = double([FLD_t0, FLD_t0(:, 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 - 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(lonGmMap + 15*pi/180*MUT - pi/2, 180/pi*colatGmMap); % Passo les matrius theta i lambda a cartesianes. x_m i y_m són matrius 91x25 (de 0º a 90º i de -180º a 180º).
    refreshdata
%    drawnow
    MAM(k) = getframe; % Agafa la imatge i la guarda en memòria.
end
movie(MAM, 1, 1)