clear all
close all

% Parametres d'entrada a modificar per l'usuari:
Field = 'jr'; % Corrents radials mesurats a una al?ada de 110 km (despr?s de mapejar les pertorbacions del camp B des de l'al?ada dels sat?l?lits (780 km) a aquesta al?ada). Jr = Jpar.r (Jpar ?s el vector corrent alineat, r ?s el vector unitari radial), per tant, positius si surten cap a fora de la ionosfera; negatius si entren.
Unitats = 'A/m^2'; % Unitats del camp Field.
Hemisf_central = 'N'; % Hemisferi central per al dibuix; 'N' o 'S'.
MaxColatProj = 60; % Graus. Maxima colatitud que es dibuixara.
PathTGCMROOT = ['..', filesep, '..', filesep];
NomJob = 'tiegcm-linux_AMPERE.job';

%
RutaiNomJob = [PathTGCMROOT, NomJob];
NomInp = [PathTGCMROOT, LlixInp(RutaiNomJob, '=', 'set input', 's', '#')];
disp(['Nom de l''arxiu d''entrada (*.inp): ', NomInp])
%NomInp = ['..', filesep, 'TIEGCM1.93', filesep, 'tiegcm1.93_default.inp'];
Start_year = LlixInp(NomInp, '=', 'START_YEAR', 'n', ';');
Start = LlixInp(NomInp, '=', 'START', 'n', ';');
Secstop = LlixInp(NomInp, '=', 'SECSTOP', 'n', ';');
Pas_temps = LlixInp(NomInp, '=', 'SECHIST', 'n', ';');
DOY_ini = Start(1);
DOY_fin = Secstop(1);
Data_ini = datestr(datenum(Start_year, 1, 1) + DOY_ini - 1, 'dd-mm-yyyy'); % Dia, mes i any en format 'DD-MM-YYYY'. Ha de correspondre al primer component de la variable START (*.inp) de l'input file del TIEGCM.
Data_fin = datestr(datenum(Start_year, 1, 1) + DOY_fin - 1, 'dd-mm-yyyy'); % Dia, mes i any en format 'DD-MM-YYYY'. Ha de correspondre al primer component de la variable START (*.inp) de l'input file del TIEGCM.
hh_ini = num2str(Start(2)); for i = length(hh_ini): 1; hh_ini = ['0', hh_ini]; end;%#ok<AGROW>
mm_ini = num2str(Start(3)); for i = length(mm_ini): 1; mm_ini = ['0', mm_ini]; end;%#ok<AGROW>
hhmm_ini = [hh_ini, ':', mm_ini];
hh_fin = num2str(Secstop(2)); for i = length(hh_fin): 1; hh_fin = ['0', hh_fin]; end;%#ok<AGROW>
mm_fin = num2str(Secstop(3)); for i = length(mm_fin): 1; mm_fin = ['0', mm_fin]; end;%#ok<AGROW>
hhmm_fin = [hh_fin, ':', mm_fin];
disp(['Data_ini = ', Data_ini, ', ', hhmm_ini])
disp(['Data_fin = ', Data_fin, ', ', hhmm_fin])
Pas_temps = Pas_temps(1)*1440 + Pas_temps(2)*60 + Pas_temps(3);
disp(['Pas_temps = ', num2str(Pas_temps), ' minuts'])
%Data = '25-08-2010'; % Dia, mes i any en format 'DD-MM-YYYY'. Ha de correspondre al primer component de la variable START (*.inp) de l'input file del TIEGCM.
%hhmm_ini = '12:00'; % Hora i minut del primer arxiu que es generar? per al TIEGCM, en format 'hh:mm'. Ha de coincidir amb els dos ?ltims components de la variable START (*.inp) de l'input file del TIEGCM. Donada la resoluci? de les dades d'AMPERE, ha de ser un m?ltiple de 120 s.
%hhmm_fin = '15:00'; % Hora i minut del darrer arxiu que es generar? per al TIEGCM, en format 'hh:mm'. Ha de coincidir amb els dos ?ltims components de la variable START (*.inp) de l'input file del TIEGCM. Donada la resoluci? de les dades d'AMPERE, ha de ser un m?ltiple de 120 s.
%Pas_temps = 20; % Minuts. Pas de temps entre arxius que es generaran per al TIEGCM. Ha de coincidir amb la variable SEWCHIST de l'input file (*.inp) del TIEGCM. Donada la resoluci? de les dades d'AMPERE, ha de ser un m?ltiple de 2 minuts.

% Obro l'arxiu netCDF, llegeixo les variables que m'interessen i tanco:
DOY_ini_ = num2str(DOY_ini); for i = length(DOY_ini_): 2; DOY_ini_ = ['0', DOY_ini_]; end; disp(['DOY_ini = ', DOY_ini_]) %#ok<AGROW>
DOY_fin_ = num2str(DOY_fin); for i = length(DOY_fin_): 2; DOY_fin_ = ['0', DOY_fin_]; end; disp(['DOY_fin = ', DOY_fin_]) %#ok<AGROW>
if DOY_ini > 1; DOY_ini_amp = DOY_ini - 1; else DOY_ini_amp = DOY_ini; disp('Compte! Falta una part del dia 001'); end
Data = datestr(datenum(Start_year, 1, 1) + DOY_ini_amp - 1, 'dd-mm-yyyy'); % Dia, mes i any en format 'DD-MM-YYYY'. Ha de correspondre al primer component de la variable START (*.inp) de l'input file del TIEGCM.
DD = Data(1: 2); MM = Data(4: 5); YYYY = Data_ini(7: 10);
RutaArxNC = ['..', filesep, 'Downloaded_data', filesep];
nomarxncN = [RutaArxNC, YYYY, MM, DD, '.north.netcdf']; % Nom de l'arxiu netCDF d'AMPERE per a l'hemisferi Nord.
nomarxncS = [RutaArxNC, YYYY, MM, DD, '.south.netcdf']; % Nom de l'arxiu netCDF d'AMPERE per a l'hemisferi Sud.
ncidN = netcdf.open(nomarxncN, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc) per a l'hemisferi Nord.
ncidS = netcdf.open(nomarxncS, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc) per a l'hemisferi Sud.
% Compte! No estic segur del significat de start_hr, start_mt i start_sc!!! ?s l'hora central de cada registre? ?s l'hora inicial? ...
UTvec =  single(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_hr'))) + single(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_mt')))/60 + single(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_sc')))/3600; % Llegeixo els arrays de temps i genero l'hora UT de cada mostra, assignant-li la variable UTvec.
UTvecS = single(netcdf.getVar(ncidS, netcdf.inqVarID(ncidS, 'start_hr'))) + single(netcdf.getVar(ncidS, netcdf.inqVarID(ncidS, 'start_mt')))/60 + single(netcdf.getVar(ncidS, netcdf.inqVarID(ncidS, 'start_sc')))/3600; % Llegeixo els arrays de temps per a l'hemisferi sud i genero l'hora UT de cada mostra, assignant-li la variable UTvecS.
if UTvec ~= UTvecS; Alto, alguna cosa no quadra; else clear UTvecS; end
FLDN = 1.0097e-6*netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, Field)); % Llegeixo la variable Jr a l'hemisferi nord i li assigno el nom FLDN. Passo a A/m2 i multiplico per un factor per passar les densitats radials de 110 km a 90 km (no caldria fer aquest darrer pas; nom?s ?s per deixar const?ncia de com caldria actuar en un cas general).
nlon = netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'nlon')); nlon = nlon(1);
nlat = netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'nlat')); nlat = nlat(1);
if nlon*nlat*length(UTvec) ~= numel(FLDN)
    Alguna cosa no quadra!
end
mlt0 = netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'mlt')); % Longitud utilitzada en el camp Field (pot ser geografica, magnetica, mlt, etc.).
mlt0(mlt0 == 24) = 0;
colatN = pi/180*netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'colat')); % Latitud N.
colatS = pi/180*netcdf.getVar(ncidS, netcdf.inqVarID(ncidS, 'colat')); % Latitud S.
netcdf.close(ncidN)
Resol_AMPERE = 60*(UTvec(2) - UTvec(1)); % Resoluci? temporal (en minuts) de les dades d'AMPERE.
if abs(Pas_temps - Resol_AMPERE*round(Pas_temps/Resol_AMPERE)) >= .01
    El pas de temps ha de ser m?ltiple de la resoluci? temporal d''AMPERE!
end
DOY_ini_amp_frac = DOY_ini_amp - 1 + UTvec(1)/24;
DOY_ini_frac = DOY_ini - 1 + (Start(2)*60 + Start(3))/1440;
DOY_fin_frac = DOY_fin - 1 + (Secstop(2)*60 + Secstop(3))/1440;
longUTvec = round((DOY_fin - DOY_ini_amp + 1)*1440/Resol_AMPERE);
UTvec = NaN(longUTvec, 1);
FLDN = NaN(size(FLDN, 1), longUTvec); FLDS = FLDN;
NumIters = round((DOY_fin_frac - DOY_ini_frac)*1440/Pas_temps + 1); if NumIters > 150; No tindras prou memoria per a la pel.licula; end
MAM(1, NumIters) = struct('cdata', [], 'colormap', []); % Inicialitzo l'estructura MAM (Movie AMpere).

Pas = 0;
for DOY = DOY_ini_amp: DOY_fin
    Pas = Pas + 1;
    Data = datestr(datenum(Start_year, 1, 1) + DOY - 1, 'dd-mm-yyyy'); % Dia, mes i any en format 'DD-MM-YYYY'. Ha de correspondre al primer component de la variable START (*.inp) de l'input file del TIEGCM.
    DD = Data(1: 2); MM = Data(4: 5);
    nomarxncN = [RutaArxNC, YYYY, MM, DD, '.north.netcdf']; % Nom de l'arxiu netCDF d'AMPERE per a l'hemisferi Nord.
    nomarxncS = [RutaArxNC, YYYY, MM, DD, '.south.netcdf']; % Nom de l'arxiu netCDF d'AMPERE per a l'hemisferi Sud.
    ncidN = netcdf.open(nomarxncN, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc) per a l'hemisferi Nord.
    ncidS = netcdf.open(nomarxncS, 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc) per a l'hemisferi Sud.
    % Compte! No estic segur del significat de start_hr, start_mt i start_sc!!! ?s l'hora central de cada registre? ?s l'hora inicial? ...
    MO = double(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_mo')));
    DY = double(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_dy')));
%    if Pas == 1; DY = DY -1; disp('Compte!!!!!!!!!!!! proves'); end
    DOY_AMP = datenum(str2double(YYYY), MO, DY) - datenum(str2double(YYYY), 1, 1) + 1;
    UTvec((Pas - 1)*1440/Resol_AMPERE + 1: Pas*1440/Resol_AMPERE)  = 24*(DOY_AMP - 1) + single(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_hr'))) + single(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_mt')))/60 + single(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_sc')))/3600; % Llegeixo els arrays de temps i genero l'hora UT de cada mostra, assignant-li la variable UTvec.
    FLDN(:, (Pas - 1)*1440/Resol_AMPERE + 1: Pas*1440/Resol_AMPERE) = 1.0097e-6*netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, Field)); % Llegeixo la variable Jr a l'hemisferi nord i li assigno el nom FLDN. Passo a A/m2 i multiplico per un factor per passar les densitats radials de 110 km a 90 km (no caldria fer aquest darrer pas; nom?s ?s per deixar const?ncia de com caldria actuar en un cas general).
    FLDS(:, (Pas - 1)*1440/Resol_AMPERE + 1: Pas*1440/Resol_AMPERE) = 1.0097e-6*netcdf.getVar(ncidS, netcdf.inqVarID(ncidS, Field)); % Llegeixo la variable Jr a l'hemisferi  sud i li assigno el nom FLDS. Passo a A/m2 i multiplico per un factor per passar les densitats radials de 110 km a 90 km (no caldria fer aquest darrer pas; nom?s ?s per deixar const?ncia de com caldria actuar en un cas general).
    netcdf.close(ncidN); netcdf.close(ncidS)
end

% M'arreglo el que he llegit a l'arxiu segons els meus interessos:
any_frac = round(10*(str2double(YYYY) + DOY_ini/365))/10;
UT_ini = (DOY_ini - 1)*24 + str2double(hhmm_ini(1: 2)) + str2double(hhmm_ini(4: 5))/60;
UT_fin = (DOY_fin - 1)*24 + str2double(hhmm_fin(1: 2)) + str2double(hhmm_fin(4: 5))/60;
PosUTvec_1 = find(abs(UTvec - UT_ini) < 0.0001);
PosUTvec_end = find(abs(UTvec - UT_fin) < 0.0001);
if isempty(PosUTvec_1) || isempty(PosUTvec_end)
    Les hores entrades no apareixen a les dades d''AMPERE!
end
PosUTvec = PosUTvec_1: round(Pas_temps/Resol_AMPERE): PosUTvec_end;
Informo = msgbox(['Es generaran ', num2str(length(PosUTvec)), ' arxius amb dades Jr. Prem Ctrl + c per abortar.'], 'Informo...', 'help'); uiwait(Informo)
PosUT = PosUTvec(1);
UT = UTvec(PosUT) - 24*floor(UTvec(PosUT)/24);
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(floor(UTvec(PosUT)/24) + 1);
FLDN_t0 = NaN(nlat, nlon); FLDS_t0 = FLDN_t0;
%FLDN_t(:) = FLDN(:, PosUTvec(1)); FLDS_t(:) = FLDS(:, PosUTvec(1));
FLDN_PosUTvec = NaN(nlat, nlon, length(PosUTvec)); FLDS_PosUTvec = FLDN_PosUTvec;
FLDN_PosUTvec(:) = FLDN(:, PosUTvec); FLDS_PosUTvec(:) = FLDS(:, PosUTvec);
FLDN_t0(:) = FLDN_PosUTvec(:, :, 1); FLDS_t0(:) = FLDS_PosUTvec(:, :, 1);

% Faig que FLDN_t i FLDS_t ocupin tot l'hemisferi corresponent.
Paslon = pi/12*(mlt0(1 + nlat) - mlt0(1)); %mltradMat(1, 2) - mltradMat(1, 1);
Paslat = colatN(2) - colatN(1);
zerosN = (max(colatN(:)) + Paslat: Paslat: pi/2)'; if zerosN(end) < pi/2; Malament; end; zerosS = (pi/2 + Paslat: Paslat: min(colatS(:)) - Paslat)'; if min(colatS(:)) - zerosS(end) > Paslat; Malament; end
FLDN_t = [FLDN_t0; zeros(length(zerosN), nlon)]; FLDS_t = [zeros(length(zerosS), nlon); FLDS_t0];
colatN = [unique(colatN); zerosN]; colatS = [zerosS; unique(colatS)];

% Defineixo la malla ?pex del TIEGCM:
[a b lonPolGm c] = Geog2Geom(any_frac, 0, 0); % ?
disp('Compte! cal posar la longitud Apex del Pol')
lonApx = pi/180*(0: 4.5: 360); %pi/12*(mltradMat - UT) - lonPolGm*pi/180;
%colatApx = pi/180*(0: 180)';
colatApx = pi/180*(90 - flipud([-90;-88.1238297810898;-86.2386381061479;-84.3344432311521;-82.4013408176051;-80.4295487348234;-78.4094762140586;-76.3318260179865;-74.1877384939343;-71.9689860630872;-69.6682254038898;-67.2793117428987;-64.7976746380798;-62.2207468739961;-59.5484273561239;-56.7835456422813;-53.9322815084405;-51.0044805887921;-48.0138007617897;-44.9776281216783;-41.9167195078328;-38.8545609882235;-35.8164740942591;-32.8285450585985;-29.9164855385726;-27.1045466399932;-24.4145973015814;-21.8654463067238;-19.4724432424165;-17.2473485865826;-15.1984257007618;-13.3306821922523;-11.6461745579382;-10.1442852711514;-8.82188355319849;-7.67329365534352;-6.69003105578526;-5.86034878972009;-5.16877182437117;-4.59593938898016;-4.11909358950423;-3.71331753226949;-3.35319031651030;-3.01415714132647;-2.67278469128132;-2.30486191304668;-1.87855705813586;-1.32758468739308;0;1.32758468739307;1.87855705813586;2.30486191304667;2.67278469128132;3.01415714132647;3.35319031651030;3.71331753226949;4.11909358950423;4.59593938898016;5.16877182437117;5.86034878972008;6.69003105578526;7.67329365534351;8.82188355319848;10.1442852711514;11.6461745579382;13.3306821922523;15.1984257007618;17.2473485865826;19.4724432424165;21.8654463067238;24.4145973015813;27.1045466399931;29.9164855385726;32.8285450585985;35.8164740942591;38.8545609882234;41.9167195078328;44.9776281216783;48.0138007617897;51.0044805887921;53.9322815084405;56.7835456422813;59.5484273561239;62.2207468739960;64.7976746380798;67.2793117428987;69.6682254038898;71.9689860630871;74.1877384939343;76.3318260179865;78.4094762140586;80.4295487348234;82.4013408176051;84.3344432311521;86.2386381061478;88.1238297810898;90]));
[lonApxMat, colatApxMat] = meshgrid(lonApx, colatApx);
lonApxMat_red = lonApxMat(:, 1: end - 1); colatApxMat_red = colatApxMat(:, 1: end - 1);
mltradMat = lonApxMat + pi/12*MLT(lonPolGm, 0, any_frac, DOY, UT, 0); % rad. No aplico la correcci? per a l'equaci? de temps (?ltima dimensi? = 0), ja que aquesta variable ser? un input del model TIEGCM (que no t? en compte la correcci?).

col1R = [0; Paslon; pi/12*min(mlt0(:)) - Paslon];
mallalon = mod(180/pi*mltradMat, 360);

% Interpolo les dades d'AMPERE a la malla Apex:
% HN:
PixelPolN = min(colatN(:)) - Paslat/2; % Dist?ncia entre el l?mit del p?xel i el pol Nord.
if PixelPolN > 0
    PixelPolN = 0;
end
col2RN = [-Paslat; 0; pi/2 - min(colatN(:)) + Paslat + PixelPolN];
mallalatN = 90 - 180/pi*(colatApxMat - PixelPolN);
RN = 180/pi*double([col1R, col2RN]);
FLDN_tApx = ltln2val(FLDN_t, RN, mallalatN, mallalon, 'bilinear');
% HS:
PixelPolS = pi - (max(colatS(:)) + Paslat/2);
if PixelPolS > 0
    PixelPolS = 0;
end
col2RS = [-Paslat; 0; pi/2 - min(colatS(:)) + Paslat - PixelPolS];
mallalatS = 90 - 180/pi*(colatApxMat + PixelPolS);
RS = 180/pi*double([col1R, col2RS]);
FLDS_tApx = ltln2val(FLDS_t, RS, mallalatS, mallalon, 'bilinear');

% Ajunto els dos hemisferis:
FLD_tApx = FLDN_tApx;
FLD_tApx(isnan(FLDN_tApx)) = FLDS_tApx(isnan(FLDN_tApx));
if isnan(FLD_tApx(1, 1))
    FLD_tApx(1, :) = mean(FLD_tApx(2, 1: end - 1))*ones(1, size(FLD_tApx, 2)); % Valor de Jr al Pol N.
end
if isnan(FLD_tApx(end, 1))
    FLD_tApx(end, :) = mean(FLD_tApx(end - 1, 1: end - 1))*ones(1, size(FLD_tApx, 2)); % Valor de Jr al Pol S.
end
if ~isempty(find(isnan(FLD_tApx), 1)); Malament; end;

% Sumo els Jr dels dos hemisferis (nord i sud) per introduir al TIEGCM de cara a solucionar per al potencial el?ctric.
numlats = round(size(FLD_tApx, 1)/2 + 0.5);
FLD_NH = FLD_tApx(1: numlats, :);
FLD_SH = FLD_tApx(numlats: end, :);
FLD_sum_NH = FLD_NH + flipud(FLD_SH);
FLD_tApx = [FLD_sum_NH; flipud(FLD_sum_NH(1: end - 1, :))]; % Poso el mateix als dos hemisferis.

% Escric el camp Jr en un arxiu tipus .\Ampere_Jr_data\Jr_YYYYDDDhhmm.dat:
RutaArxJr = ['..', filesep, 'Ampere_Jr_data', filesep];
NomArxJr = [RutaArxJr, 'Jr_', YYYY, DOY_, hh, mm, '.dat'];
Arx_Jr = fopen(NomArxJr, 'w'); % Obro l'arxiu per a escriptura
fprintf(Arx_Jr, [repmat('%+1.3e ', 1, length(colatApx)), '\n'], flipud(FLD_tApx)); % Compte! L'arxiu generat cont? una matriu (nmlon + 1) x nmlat (81x97), on les files s?n longituds i les columnes s?n latituds. Latitud -90 ?s a la columna 1, ..., latitud 90 ?s a l'?ltima columna.
fclose(Arx_Jr);

% Estableixo els nivells de colors adequats i els seus limits de cara a dibuixar:
FLD2 = FLDN_PosUTvec + flipdim(FLDS_PosUTvec, 1); % Compte! No est? ben fet. Nom?s afecta els l?mits de colors del dibuix.
minFLD = min(FLD2(:)); maxFLD = max(FLD2(:));
maxFLD = max(abs(maxFLD), abs(minFLD)); minFLD = -maxFLD;
%minFLD=-1e-6; maxFLD=1e-6;
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];

% Creo un mapa polar:
figure('Position', [300, 50, 800, 600])
h = polar(-pi: pi/100: pi, MaxColatProj/2*ones(1, 201), '-k'); delete(h)
title({['J_r component of FACs (positive upward) - ', Hemisf_central, 'H projection']; 'Jr_N_H + Jr_S_H'}, 'FontSize', 13)
xlim([-MaxColatProj MaxColatProj])
ylim([-MaxColatProj MaxColatProj])
xlabel('00 MLT', 'FontSize', 12);
ylabel('18 MLT', 'FontSize', 12);
hb = colorbar('FontSize', 12);
ylabel(hb, Unitats, 'FontSize', 12, 'Rotation', 0, 'Units', 'centimeters', 'Position', [2.28, 13.25, 0])
hold on

% Afegeixo l?nies de costa. Compte! El mapa est? en un sistema de coordenades geomagn?tiques; no Apex. Aix? nom?s afecta el dibuix.
load coast
[lonGmMap latGmMap lonPolGm colatPolGm] = Geog2Geom(any_frac, long, lat); % ?. Converteixo la malla geogr?fica a geomagn?tica.
clear long lat
lonGmMap = lonGmMap*pi/180; %rad
colatGmMap = pi/180*(90 - latGmMap); %rad
MUT = MLT(lonPolGm, 0, any_frac, DOY, UT, 0); % No aplico la correcci? per a l'equaci? de temps (?ltima dimensi? = 0), ja que aquesta variable ser? un input del model TIEGCM (que no t? en compte la correcci?).
if Hemisf_central == 'N' % HN
    colatApxMat(colatApxMat > pi/180*MaxColatProj) = NaN;
    [x_0, y_0] = pol2cart(mltradMat - pi/2, 180/pi*colatApxMat); % Passo les matrius de colatitud i MLT a cartesianes.
    colatGmMap(colatGmMap > pi/180*MaxColatProj) = NaN;
    [x_m, y_m] = pol2cart(lonGmMap + pi/12*MUT - pi/2, 180/pi*colatGmMap); % Passo les coordenades de les costes a cartesianes.
else % HS
    colatApxMat(colatApxMat < pi - pi/180*MaxColatProj) = NaN;
    [x_0, y_0] = pol2cart(mltradMat - pi/2, 180 - 180/pi*colatApxMat); % Passo les matrius de colatitud i MLT a cartesianes.
    colatGmMap(colatGmMap < pi - pi/180*MaxColatProj) = NaN;
    [x_m, y_m] = pol2cart(lonGmMap + pi/12*MUT - pi/2, 180 - 180/pi*colatGmMap); % Passo les coordenades de les costes a cartesianes.
end

% Dibuixo les dades de Jr en un mapa polar:
[C, hc] = contourf(x_0, y_0, FLD_tApx, '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_tApx')

hmap = plot(x_m, y_m, 'b'); % Dibuixo les l?nies de costa.
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(floor(UTvec(PosUT)/24) + datenum(str2double(YYYY), 1, 1)); ['DOY = ' DOY_]}, 'linestyle', 'none', 'FontSize', 12);

% Itero per a cada pas de temps:
k = 1;
MAM(k) = getframe; % Agafa la imatge i la guarda en mem?ria.
for PosUT = PosUTvec(2: end)
    k = k + 1;
    UT = UTvec(PosUT) - 24*floor(UTvec(PosUT)/24);
    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>
    DOY_ = num2str(floor(UTvec(PosUT)/24) + 1);
    MUT = MLT(lonPolGm, 0, any_frac, DOY, UT, 0);
    FLDN_t0 = NaN(nlat, nlon); FLDS_t0 = FLDN_t0;
    FLDN_t0(:) = FLDN_PosUTvec(:, :, k); FLDS_t0(:) = FLDS_PosUTvec(:, :, k);
    FLDN_t = [FLDN_t0; zeros(length(zerosN), nlon)]; FLDS_t = [zeros(length(zerosS), nlon); FLDS_t0];
    mltradMat = lonApxMat + pi/12*MLT(lonPolGm, 0, any_frac, DOY, UT, 0); % rad. No aplico la correcci? per a l'equaci? de temps (?ltima dimensi? = 0), ja que aquesta variable ser? un input del model TIEGCM (que no t? en compte la correcci?).
    mallalon = mod(180/pi*mltradMat, 360);
    FLDN_tApx = ltln2val(FLDN_t, RN, mallalatN, mallalon, 'bilinear'); % Interpolo HN.
    FLDS_tApx = ltln2val(FLDS_t, RS, mallalatS, mallalon, 'bilinear'); % Interpolo HS.
    
    % Ajunto els dos hemisferis:
    FLD_tApx = FLDN_tApx;
    FLD_tApx(isnan(FLDN_tApx)) = FLDS_tApx(isnan(FLDN_tApx));
    if isnan(FLD_tApx(1, 1))
        FLD_tApx(1, :) = mean(FLD_tApx(2, 1: end - 1))*ones(1, size(FLD_tApx, 2)); % Valor de Jr al Pol N.
    end
    if isnan(FLD_tApx(end, 1))
        FLD_tApx(end, :) = mean(FLD_tApx(end - 1, 1: end - 1))*ones(1, size(FLD_tApx, 2)); % Valor de Jr al Pol S.
    end
    if ~isempty(find(isnan(FLD_tApx), 1)); Malament; end;
    
    % Sumo els Jr dels dos hemisferis (nord i sud) per introduir al TIEGCM de cara a solucionar per al potencial el?ctric.
    FLD_NH = FLD_tApx(1: numlats, :);
    FLD_SH = FLD_tApx(numlats: end, :);
    FLD_sum_NH = FLD_NH + flipud(FLD_SH);
    FLD_tApx = [FLD_sum_NH; flipud(FLD_sum_NH(1: end - 1, :))]; % Poso el mateix als dos hemisferis.
    
    % Escric el camp Jr en un arxiu tipus .\Ampere_Jr_data\Jr_YYYYDDDhhmm.dat:
    NomArxJr = [RutaArxJr, 'Jr_', YYYY, DOY_, hh, mm, '.dat'];
    Arx_Jr = fopen(NomArxJr, 'w'); % Obro l'arxiu per a escriptura
    fprintf(Arx_Jr, [repmat('%+1.3e ', 1, length(colatApx)), '\n'], flipud(FLD_tApx)); % Compte! L'arxiu generat cont? una matriu (nmlon + 1) x nmlat (81x97), on les files s?n longituds i les columnes s?n latituds. Latitud -90 ?s a la columna 1, ..., latitud 90 ?s a l'?ltima columna.
    fclose(Arx_Jr);
    
    % Passo les matrius de colatitud i MLT a cartesianes:
    if Hemisf_central == 'N' % HN
        [x_0, y_0] = pol2cart(mltradMat - pi/2, 180/pi*colatApxMat); % Passo les matrius de colatitud i MLT a cartesianes.
        [x_m, y_m] = pol2cart(lonGmMap + pi/12*MUT - pi/2, 180/pi*colatGmMap); % Passo les coordenades de les costes a cartesianes.
    else % HS
        [x_0, y_0] = pol2cart(mltradMat - pi/2, 180 - 180/pi*colatApxMat); % Passo les matrius de colatitud i MLT a cartesianes.
        [x_m, y_m] = pol2cart(lonGmMap + pi/12*MUT - pi/2, 180 - 180/pi*colatGmMap); % Passo les coordenades de les costes a cartesianes.
    end
    
    delete(han)
    han = annotation('textbox', [0.18 0.11 .1 .1], 'String', {[hh, ':', mm, ' UT']; datestr(floor(UTvec(PosUT)/24) + datenum(str2double(YYYY), 1, 1)); ['DOY = ' DOY_]}, 'linestyle', 'none', 'FontSize', 12);
    refreshdata
    MAM(k) = getframe; % Agafa la imatge i la guarda en mem?ria.
end
movie(MAM, 1, 1)