%
% 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:
%Posilev = 29; % Posici? dintre del vector ilev (d?na l'al?ada com un nivell de pressi? constant dintre de la ionosfera).
PosUTvec = 1: 3; % Posici? dintre del vector UT per al que es dibuixar? el potencial.
Field = 'JQR';
Unitats = 'A/m^2';
PathTGCMROOT = ['..', filesep, '..', filesep];
LastOrOther = 'last'; % Cal posar 'last' o 'other', depenent de si es vol obrir l'?ltim arxiu SECOUT generat pel TIEGCM (el buscar? a la carpeta RutaArxNC_last), o b? un altre; per ex., un guardat a la carpeta TIEGCM_results.
% En cas que s'hagi triat 'last':
RutaArxNC_last = [PathTGCMROOT, 'tiegcm-linux', filesep]; % Path on el TIEGCM genera els arxius SECOUT.
%NomJob = 'tiegcm-linux_AMPERE.job'; % Nom de l'arxiu .job on buscar l'arxiu .inp que contindr? el nom de l'arxiu SECOUT en cas que s'hagi triat l'opci? 'last'. No t? cap efecte en cas que s'hagi triat 'other'.
NomJob = 'tiegcm-linux_AMPERE.job'; % Nom de l'arxiu .job on buscar l'arxiu .inp que contindr? el nom de l'arxiu SECOUT en cas que s'hagi triat l'opci? 'last'. No t? cap efecte en cas que s'hagi triat 'other'.
% En cas que s'hagi triat 'other':
RutaiNomArxNC = ['..', filesep, 'TIEGCM_results', filesep, '*.nc']; % Ruta i nom de l'arxiu SECOUT en cas que s'hagi triat l'opci? 'other'. No t? cap efecte en cas que s'hagi triat 'last'.
MaxColatProj = 60; % ?. M?nima colatitud de la projecci?.
%dimlon = 'lon'; % Dimensi? longitud utilitzada pel camp Field. Cal fer un 'ncdump -c nomarxnc' (a Linux) per veure quina dimensi? utilitza la variable Field!
%dimlon = 'mlon'; % 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!
%dimlat = 'mlat'; % Dimensi? latitud utilitzada pel camp Field. Cal fer un 'ncdump -c nomarxnc' (a Linux) per veure quina dimensi? utilitza la variable Field!

trobat_altura = false;
if strcmp(LastOrOther, 'last')
    RutaiNomJob = [PathTGCMROOT, NomJob];
    RutaiNomInp = [PathTGCMROOT, LlixInp(RutaiNomJob, '=', 'set input', 's', '#')];
    NomArxSECOUT = LlixInp(RutaiNomInp, '=', 'SECOUT', 's', ';');
    NomArxSECOUT = NomArxSECOUT(2: end - 1);
    if strfind(NomArxSECOUT, ''''); SECOUT ambigu; end
    RutaiNomArxNC =  [RutaArxNC_last, NomArxSECOUT];
end
disp(['Ruta i nom de l''arxiu SECOUT (*.nc) a dibuixar: ', RutaiNomArxNC])

% Obro l'arxiu netCDF, llegeixo les variables que m'interessen i tanco:
ncid = netcdf.open(RutaiNomArxNC, '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 de la variable ', Field, ':'])
for dimid = dimids
    [dimnom, dimllarg] = netcdf.inqDim(ncid, dimid);
    disp([dimnom, ', ', num2str(dimllarg)])
    if strcmp(dimnom, 'ilev');
        trobat_altura = true;
        prompt = {'Enter height index to be drawn:'};
        dlg_title = 'Input for drawing height';
        num_lines = 1;
        def = {'1'};
        Posilev = str2double(cell2mat(inputdlg(prompt,dlg_title,num_lines,def)));
        ilev = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ilev')); ilev = ilev(Posilev); % Llegeixo el vector ilev.
    end
end
if length(dimids) == 4 && ~trobat_altura; Les dimensions no em quadren; end
any = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'year')); any = double(any(1));
DOY = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'mtime')); DOY = double(DOY(1));
any_frac = any + DOY/365;
UTvec = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ut')); % Llegeixo el vector ut i li assigno la variable UTvec.
FLD = netcdf.getVar(ncid, netcdf.inqVarID(ncid, Field)); % Llegeixo la variable especificada i li assigno la variable FLD.
FLD(FLD > 1e36) = NaN;
%lon_coor = pi/180*netcdf.getVar(ncid, netcdf.inqVarID(ncid, dimlon))'; % Longitud utilitzada en el camp Field (pot ser geografica, magnetica, mlt, etc.).
lon_coor = pi/180*netcdf.getVar(ncid, dimids(1))'; % Podria automatitzar-ho aix?, per? m'estimo m?s assegurar-me fent un ncdump i veient de quines dimensions depen la variable Field.
%colat_coor = pi/180*(90 - wrev(netcdf.getVar(ncid, netcdf.inqVarID(ncid, dimlat))));  % Longitud utilitzada en el camp Field (pot ser geografica, magnetica, mlt, etc.). Cal invertir el vector de latituds pr?viament, ja que comen?a a l'hemisferi Sud.
colat_coor = pi/180*(90 - wrev(netcdf.getVar(ncid, dimids(2))));  % Longitud utilitzada en el camp Field (pot ser geografica, magnetica, mlt, etc.). Cal invertir el vector de latituds pr?viament, ja que comen?a a l'hemisferi Sud.
netcdf.close(ncid)

% M'arreglo el que he llegit a l'arxiu segons els meus interessos:
PosUT = PosUTvec(1);
UT = UTvec(PosUT); % UT ?s ara el temps universal pel qual es far? el dibuix.
hh = num2str(floor(UT)); mm = num2str(round(60*(UT - floor(UT))));
if length(hh) == 1; hh = ['0', hh]; end; if length(mm) == 1; mm = ['0', mm]; end
DOY_ = num2str(DOY);
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 = colat_coor(1: floor((length(colat_coor) + 1)/2)); % M'interessa nom?s l'hemisferi Nord.
colat_coor(colat_coor > pi/180*MaxColatProj) = NaN;
[lonMat, colatMat] = meshgrid(lon_coor, colat_coor); % Genero la malla.
FLD_t_ilev0 = NaN(size(FLD, 1), size(FLD, 2));
FLD_PosUTvec = NaN(size(FLD, 1), size(FLD, 2), length(PosUTvec));
if ~trobat_altura
    FLD_PosUTvec(:) = FLD(:, :, PosUTvec); % M'interessa el potencial per a una al?ada determinada (donada per ilev) i per a un temps UT determinat.
else
    FLD_PosUTvec(:) = FLD(:, :, Posilev, PosUTvec); % M'interessa el potencial per a una al?ada determinada (donada per ilev) i per a un temps UT determinat.
end
FLD_t_ilev0(:) = FLD_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: floor((size(FLD_t_ilev1, 1) + 1)/2), :); % M'interessa nom?s l'hemisferi Nord.
FLD_t_ilev2 = FLD_t_ilev1;
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.

% Estableixo els nivells de colors adequats i els seus limits:
minFLD = min(FLD_PosUTvec(:)); maxFLD = max(FLD_PosUTvec(:));
maxFLD = max(abs(maxFLD), abs(minFLD)); minFLD = -maxFLD;
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 = [-5, 30];
%Salt = 1;
%LlistaNivells = LimitsFLD(1): Salt: LimitsFLD(2);

% Creo un mapa polar:
figure('Position', [300, 90, 800, 600])
h = polar(-pi: pi/100: pi, MaxColatProj/2*ones(1, 201), '-k'); delete(h)
title({Field}, 'FontSize', 13)
xlim([-MaxColatProj MaxColatProj])
ylim([-MaxColatProj MaxColatProj])
if dimids(1) == 1 % Coordenades geografiques.
    xlabel('00 LT', 'FontSize', 12);
    ylabel('18 LT', 'FontSize', 12);
elseif dimids(1) == 5 % Coordenades geomagnetiques dipolars.
    xlabel('00 MLT', 'FontSize', 12);
    ylabel('18 MLT', 'FontSize', 12);
end
hb  = colorbar('FontSize', 12);
ylabel(hb, Unitats, 'FontSize', 12, 'Rotation', 0, 'Units', 'centimeters', 'Position', [2.3, 13.6, 0])
hold on

% Passo les matrius theta i lambda a cartesianes:
if dimids(1) == 1 % Coordenades geografiques.
    [x_0, y_0] = pol2cart(lonMat + pi/12*UT - pi/2, 180/pi*colatMat); % x_0 i y_0 s?n matrius 91x25 (de 0? a 90? i de -180? a 180?).    
elseif dimids(1) == 5 % Coordenades geomagnetiques dipolars.
    [a b lonPolGm c] = Geog2Geom(any_frac, 0, 0); % ?. Converteixo la malla geogr?fica a geomagn?tica.
    MUT = MLT(lonPolGm, 0, any_frac, DOY, UT, 1);
    [x_0, y_0] = pol2cart(lonMat + pi/12*MUT - pi/2, 180/pi*colatMat); % x_0 i y_0 s?n matrius 91x25 (de 0? a 90? i de -180? a 180?).
else
    Error, no se quina longitud utilitzes
end

% Dibuixo les dades de FLD en un mapa polar:
[C, hc] = contourf(x_0, y_0, FLD_t_ilev, 'LineStyle', 'none', 'DisplayName', [Field, '(', 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_ilev')

% Afegeixo l?nies de costa.
load coast
if dimids(1) == 1 % Coordenades geografiques.
    lonMap = pi/180*long; % rad
    colatMap = pi/180*(90 - lat); % rad
    colatMap(colatMap > pi/180*MaxColatProj) = NaN;
    [x_m, y_m] = pol2cart(lonMap + pi/12*UT - pi/2, 180/pi*colatMap);
elseif dimids(1) == 5 % Coordenades geomagnetiques dipolars.
    [lonMap latMap lonPolGm colatPolGm] = Geog2Geom(any_frac, long, lat); % ?. Converteixo la malla geogr?fica a geomagn?tica.
    lonMap = lonMap*pi/180; %rad
    colatMap = pi/180*(90 - latMap); %rad
    colatMap(colatMap > pi/180*MaxColatProj) = NaN;
    MUT = MLT(lonPolGm, 0, any_frac, DOY, UT, 1);
    [x_m, y_m] = pol2cart(lonMap + pi/12*MUT - pi/2, 180/pi*colatMap);
end
clear long lat
hmap = plot(x_m, y_m, 'k');
set(hmap, 'XDataSource', 'x_m', 'YDataSource', 'y_m')

% Pinto meridians i paral.lels:
for i = 0: 30: 359
    polar([i i]*pi/180, [0 MaxColatProj],  ':k') % Meridians
end
for i = 20: 20: MaxColatProj
    polar(-pi: pi/100: pi, i*ones(1, 201), ':k') % Paral.lels
end
han =  annotation('TextBox', [0.18 0.11 .1 .1], 'String', {[hh, ':', mm, ' UT']; datestr(DOY - 1 + datenum(any, 1, 1)); ['DOY = ' DOY_]}, 'LineStyle', 'none', 'FontSize', 12);
if trobat_altura
    han2 = annotation('TextBox', [.65 0.11 .1 .1], 'String', {[num2str(90 - MaxColatProj), '^o lat.']; ['ln(p0/p) = ', num2str(ilev)]}, 'LineStyle', 'none', 'FontSize', 12);
else
    han2 = annotation('TextBox', [.65 0.11 .1 .1], 'String', {[num2str(90 - MaxColatProj), '^o lat.']}, 'LineStyle', 'none', 'FontSize', 12);    
end

% Itero per a cada pas de temps:
MT(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.
    hh = num2str(floor(UT)); mm = num2str(round(60*(UT - floor(UT))));
    if length(hh) == 1; hh = ['0', hh]; end; if length(mm) == 1; mm = ['0', mm]; end %#ok<AGROW>
    FLD_t_ilev0(:) = FLD_PosUTvec(:, :, k); % 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: floor((size(FLD_t_ilev1, 1) + 1)/2), :); % M'interessa nom?s l'hemisferi Nord.
    FLD_t_ilev2 = FLD_t_ilev1;
    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.
    
    if dimids(1) == 1 % Coordenades geografiques.
        [x_0, y_0] = pol2cart(lonMat + pi/12*UT - pi/2, 180/pi*colatMat); % x_0 i y_0 s?n matrius 91x25 (de 0? a 90? i de -180? a 180?).    
        [x_m, y_m] = pol2cart(lonMap + pi/12*UT - pi/2, 180/pi*colatMap);     
    elseif dimids(1) == 5 % Coordenades geomagnetiques dipolars.
        MUT = MLT(lonPolGm, 0, any_frac, DOY, UT, 1); % Aplico la correcci? per a l'equaci? de temps (?ltima dimensi? = 1).
        [x_0, y_0] = pol2cart(lonMat + pi/12*MUT - 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 + pi/12*MUT - pi/2, 180/pi*colatMap); % 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?).
    end
    
    delete(han)
    han = annotation('TextBox', [0.18 0.11 .1 .1], 'String', {[hh, ':', mm, ' UT']; datestr(DOY - 1 + datenum(any, 1, 1)); ['DOY = ' DOY_]}, 'LineStyle', 'none', 'FontSize', 12);
    refreshdata
    MT(k) = getframe; % Agafa la imatge i la guarda en mem?ria.
end
movie(MT, 3, 1)