clear all
close all

% Paràmetres 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.
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 = '24:08'; % 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 = '24:08'; % 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 = 2; % 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.
Hemisf_central = 'N'; % Hemisferi central per al dibuix; 'N' o 'S'.
MaxColatProj = 60; % º. Màxima colatitud que es dibuixarà.

% Obro l'arxiu netCDF, llegeixo les variables que m'interessen i tanco:
DD = Data(1: 2); MM = Data(4: 5); YYYY = Data(7: 10);
DOY = daysact(['01-01-', YYYY], [MM, '-', DD, '-', YYYY]) + 1; % Ha de ser el primer component de la variable START de l'input file (*.inp) del TIEGCM.
DOY_ = num2str(DOY); for i = length(DOY_): 2; DOY_ = ['0', DOY_]; end; disp(['DOY = ', DOY_]) %#ok<AGROW>
nomarxncN = ['.', filesep, YYYY, MM, DD, '.north.nc']; % Nom de l'arxiu netCDF d'AMPERE per a l'hemisferi Nord.
nomarxncS = ['.', filesep, YYYY, MM, DD, '.south.nc']; % 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 = 24*(single(netcdf.getVar(ncidN, netcdf.inqVarID(ncidN, 'start_dy'))) - str2double(DD)) + 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 = 24*(single(netcdf.getVar(ncidS, netcdf.inqVarID(ncidN, 'start_dy'))) - str2double(DD)) + 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
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
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).
FLDS = 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).
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); netcdf.close(ncidS)

% M'arreglo el que he llegit a l'arxiu segons els meus interessos:
any_frac = round(10*(str2double(YYYY) + DOY/365))/10;
UT_ini = str2double(hhmm_ini(1: 2)) + str2double(hhmm_ini(4: 5))/60;
UT_fin = 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)        
UT = UTvec(PosUTvec(1));
FLDN_t = NaN(nlat, nlon); FLDS_t = FLDN_t;
%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_t(:) = FLDN_PosUTvec(:, :, 1); FLDS_t(:) = 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_t; zeros(length(zerosN), nlon)]; FLDS_t = [zeros(length(zerosS), nlon); FLDS_t];
colatN = [unique(colatN); zerosN]; colatS = [zerosS; unique(colatS)];

% Divideixo per -sin(I) per a passar corrents radials a FACs (positius en el sentit del camp magnètic). Però sembla que no cal, ja que segons Astrid, l'entrada del TIEGCM són corrents radials.
%sinI_N = repmat(2*cos(colatN)./sqrt(1 + 3*cos(colatN).^2), 1, nlon); sinI_S = repmat(2*cos(colatS)./sqrt(1 + 3*cos(colatS).^2), 1, nlon);
%FLDN_t = -FLDN_t./(sinI_N + eps); FLDS_t = -FLDS_t./(sinI_S + eps); % Desfaig la projecció radial. Ara tinc FACs.

%Filtro amb una mitjana mòbil al llarg de la longitud d'ona més curta que contenen els harmònics. Ho faig al llarg de les latituds:
nmax = 60;
mmax = 5;
%nemmlat = 2*round((ceil(2*pi/nmax/Paslat) - 1)/2) - 1; % Número d'elements de la mitjana mòbil. Depèn de nmax, ja que 360/nmax és la longitud d'ona més curta que em trobaré.
nemmlat = 1; % No filtro.
filtlat = ones(nemmlat, 1); % Filtre.
ExisteixN = FLDN_t; ExisteixN(:) = 1; ExisteixS = FLDS_t; ExisteixS(:) = 1;
FLDN_t_filt = filter2(filtlat, FLDN_t)./filter2(filtlat, ExisteixN); FLDS_t_filt = filter2(filtlat, FLDS_t)./filter2(filtlat, ExisteixS);
FLDN_t = FLDN_t_filt; FLDS_t = FLDS_t_filt;

% 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 Àpex:
% 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.

% Filtro!!! Compte! No puc filtrar la malla Apex de TIEGCM, ja que no és regular.
%NumeradorFiltre = [0.000837164591190549,0.000630088428244303,0.000117016843802876,-0.00121993975387921,-0.00332542341352173,-0.00574692952437155,-0.00769394282840493,-0.00829077778658227,-0.00695524134107235,-0.00374942375230463,0.000481574780048021,0.00430483829567512,0.00621337352781100,0.00528242170471528,0.00168685617863516,-0.00319545745934426,-0.00722468309056643,-0.00836788038727407,-0.00565972803429091,0.000179656876472077,0.00682755832844448,0.0112162976865154,0.0108705945822515,0.00519156780704766,-0.00396030394779504,-0.0127427426864552,-0.0168337075123898,-0.0134048554494479,-0.00272688786944293,0.0113997298499187,0.0227782008694017,0.0252543616580168,0.0155816632292830,-0.00447728139280901,-0.0278377722520717,-0.0439625206758760,-0.0424408414544781,-0.0170528693128577,0.0312321802089387,0.0935906806841477,0.155466371458946,0.200894019011274,0.217573425599882,0.200894019011274,0.155466371458946,0.0935906806841477,0.0312321802089387,-0.0170528693128577,-0.0424408414544781,-0.0439625206758760,-0.0278377722520717,-0.00447728139280901,0.0155816632292830,0.0252543616580168,0.0227782008694017,0.0113997298499187,-0.00272688786944293,-0.0134048554494479,-0.0168337075123898,-0.0127427426864552,-0.00396030394779504,0.00519156780704766,0.0108705945822515,0.0112162976865154,0.00682755832844448,0.000179656876472077,-0.00565972803429091,-0.00836788038727407,-0.00722468309056643,-0.00319545745934426,0.00168685617863516,0.00528242170471528,0.00621337352781100,0.00430483829567512,0.000481574780048021,-0.00374942375230463,-0.00695524134107235,-0.00829077778658227,-0.00769394282840493,-0.00574692952437155,-0.00332542341352173,-0.00121993975387921,0.000117016843802876,0.000630088428244303,0.000837164591190549;];
%NumeradorFiltre = [0.00181823971512725,0.00517005788790013,0.00856960842158199,0.0103199378817500,0.00757480877420618,0.00102502966908523,-0.00651062351066224,-0.00961782382684962,-0.00551247556294510,0.00421798998168117,0.0125298811090595,0.0122627325941355,0.00131221813989646,-0.0135787089102958,-0.0207727129892657,-0.0118248304707069,0.0100426666260028,0.0297437433025841,0.0292807613287389,0.00218309677238062,-0.0376784936076520,-0.0600760175020828,-0.0358038956666434,0.0431439617305824,0.153584726816896,0.250030168228445,0.288240429044952,0.250030168228445,0.153584726816896,0.0431439617305824,-0.0358038956666434,-0.0600760175020828,-0.0376784936076520,0.00218309677238062,0.0292807613287389,0.0297437433025841,0.0100426666260028,-0.0118248304707069,-0.0207727129892657,-0.0135787089102958,0.00131221813989646,0.0122627325941355,0.0125298811090595,0.00421798998168117,-0.00551247556294510,-0.00961782382684962,-0.00651062351066224,0.00102502966908523,0.00757480877420618,0.0103199378817500,0.00856960842158199,0.00517005788790013,0.00181823971512725;];
%FLD_tApx = filter2(NumeradorFiltre', FLD_tApx);

%FLD_tApx_red = FLD_tApx(:, 1: end - 1);

% Escric el camp Jr en un arxiu tipus .\ampere_data\Jr_YYYYDDDhhmm.dat:
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
nomarx = ['.', filesep, 'ampere_data', filesep, 'Jr_', YYYY, DOY_, hh, mm, '.dat'];
Arx_Jr = fopen(nomarx, '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);

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

% Estableixo els nivells de colors adequats i els seus limits de cara a dibuixar:
%FLDN2 = FLDN(:, PosUTvec); FLDS2 = FLDS(:, PosUTvec);
%minFLDN = min(FLDN2(:)); minFLDS = min(FLDS2(:)); maxFLDN = max(FLDN2(:)); maxFLDS = max(FLDS2(:));
%minFLD = min(minFLDN, minFLDS); maxFLD = max(maxFLDN, maxFLDS);
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(datenum([MM, '-', DD, '-', YYYY], 'mm-dd-yyyy')); ['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); % UT és ara el temps universal pel qual es farà el dibuix.
    MUT = MLT(lonPolGm, 0, any_frac, DOY, UT, 0);
    FLDN_t0 = NaN(nlat, nlon); FLDS_t0 = FLDN_t0;
%    FLDN_t0(:) = FLDN(:, PosUT); FLDS_t0(:) = FLDS(:, PosUT);
    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_data\Jr_YYYYDDDhhmm.dat:
    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>
    nomarx = ['.', filesep, 'ampere_data', filesep, 'Jr_', YYYY, DOY_, hh, mm, '.dat'];
    Arx_Jr = fopen(nomarx, '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(datenum([MM, '-', DD, '-', YYYY], 'mm-dd-yyyy')); ['DOY = ' DOY_]}, 'linestyle', 'none', 'FontSize', 12);
    refreshdata
    MAM(k) = getframe; % Agafa la imatge i la guarda en memòria.
end
movie(MAM, 1, 1)