function MLT = MLT(LonObs0, LatObs0, any, DOY, UT, Corr_Eq_t) % Càlcul de MLT en funció de les coordenades geogràfiques del punt desitjat, any (fraccionari), dia de l'any i UT. Compte! es basa en coordenades geomagnètiques; no magnètiques.
% UT és el temps solar mitjà de Greenwich. Es corregeix per l'equació de temps per a trobar el punt subsolar.
% La Declinació solar es corregeix per una òrbita el·líptica.
LonObs = LonObs0; LonObs = LonObs(:); MLT = LonObs0;
LatObs = LatObs0; LatObs = LatObs(:);

%colatPolGm = 90 - 79.744; % època 2005.0
%lonPolGm = -71.784;
RutaiNomArxIGRF = './IGRF11.txt';
disp(['Utilitzant ', RutaiNomArxIGRF])
% Import the file
newData1 = importdata(RutaiNomArxIGRF);
data = newData1.data;
textdata = newData1.textdata;
clear newData1
CoefsIGRF = data(:, 3: end - 1);
SV = data(:, end);
CapFil = textdata(1, 4: end);
AnysCapFil = str2num(cell2mat(CapFil(1: end-1)'))';
prany = str2double(CapFil{1});
ultany = str2double(CapFil{end - 1});
CapCol = [textdata(2: end, 1), num2cell(data(:, 1: 2))];
if (any >= prany) && (any <= ultany)
    g10 = interp1(AnysCapFil', CoefsIGRF(1, :)', any);
    g11 = interp1(AnysCapFil', CoefsIGRF(2, :)', any);
    h11 = interp1(AnysCapFil', CoefsIGRF(3, :)', any);
elseif any < prany
    disp('Cal refinar el model IGRF!')
    g10 = interp1(AnysCapFil', CoefsIGRF(1, :)', any, 'linear', 'extrap');
    g11 = interp1(AnysCapFil', CoefsIGRF(2, :)', any, 'linear', 'extrap');
    h11 = interp1(AnysCapFil', CoefsIGRF(3, :)', any, 'linear', 'extrap');
elseif any > ultany
    if any > ultany + 5
        disp('Compte, cal actualitzar el model IGRF!')
    end
    g10 = CoefsIGRF(1, end) + SV(1)*(any - ultany);
    g11 = CoefsIGRF(2, end) + SV(2)*(any - ultany);
    h11 = CoefsIGRF(3, end) + SV(3)*(any - ultany);
end
colatPolGm = acosd(-g10/sqrt(g10^2 + g11^2 + h11^2));
lonPolGm = 180/pi*atan2(-h11, -g11);
disp(['Longitud geogràfica del pol Nord geomagnètic (dipol): ', num2str(lonPolGm), 'º'])

DecSol = asind(sind(23.45)*sind(360/365.25*(DOY-80))); % Declinació solar. Òrbita circular a velocitat constant per l'eclíptica (23.45º).
DecSol = DecSol - (4.182*sin(0.009179*DOY + 4.317) + 4.41*sin(0.01094*DOY + 0.9783) + 0.3535*sin(0.03456*DOY + 1.743)); % Correcció òrbita el·líptica.
% Equació del temps (temps solar mitjà - temps solar vertader) aproximada (en hores!). Aquesta correcció obeeix a 1- Un angle de l'eclíptica (23.45º) diferent de zero, i 2- a l'el·lipticitat de l'òrbita terrestre:
if Corr_Eq_t
    ET = 24/(2*pi)*(atan2(sind(360/365.25*(DOY - 80)), cosd(360/365.25*(DOY - 80))) - atan2(sind(360/365.25*(DOY - 80)), cosd(23.45)*cosd(360/365.25*(DOY - 80)))) + 8/60*sind(360/365.25*(DOY-3));
    UT = UT - ET; % Ara tinc el temps solar vertader (o aparent) de Greenwich.
end
%LatObs = 45;
%LonObs = 0;
%UT = 0;
%UT=-lonPolGm/15+6;
%colatPolGm=0;
% lonPolGm=-72;
% colatPolGm=10;
% DecSol = 0;
%LonObs = lonPolGm;
[XPol, YPol, ZPol] = sph2cart(pi/180*(lonPolGm + 15*UT), pi/2 - pi/180*colatPolGm, 1);
VecPol = [XPol, YPol, ZPol]';
[XSol, YSol, ZSol] = sph2cart(pi, pi/180*DecSol, 1);
VecSol = [XSol, YSol, ZSol]';
[XObs, YObs, ZObs] = sph2cart(pi/180*(LonObs + 15*UT), pi/180*LatObs, 1);
VecObs =[XObs, YObs, ZObs]';
VecPolSol = cross(VecPol, -VecSol); VecPolSol = VecPolSol/sqrt(VecPolSol(1)^2 + VecPolSol(2)^2 + VecPolSol(3)^2);
VecPol = repmat(VecPol, 1, length(LonObs));
VecPolObs = cross(VecPol, VecObs); VecPolObs = VecPolObs./repmat(sqrt(sum(VecPolObs.^2)), 3, 1);
MLT0 = acos(VecPolSol'* VecPolObs)*12/pi;
VecPolSol = repmat(VecPolSol, 1, length(LonObs));
PV = cross(VecPolSol, VecPolObs);
PVZ = sign(PV(3, :));
MLT0 = mod(MLT0.*PVZ, 24);
MLT(:) = MLT0;