% Dibuixo el potencial electrost?tic a partir d'un arxiu netCDF.
clear all
close all
PosUT = 30; % Posici? dintre del vector ut.
Posilev = 15; % Posici? dintre del vector ilev (d?na l'al?ada com un nivell de pressi? constant dintre de la ionosfera).
ncid = netcdf.open('/hao/aim2/marsal/GangLu/amie/sep16_2005_nh.nc', 'NC_NOWRITE'); % Obro l'arxiu netCDF (= *.nc)
UT = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'ut')); % Llegeixo el vector ut i li assigno la variable UT
UT = UT(PosUT); % UT ?s ara el temps universal pel qual es far? el dibuix.
Phi = netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'POTEN')); % Llegeixo la matriu POTEN i li assigno la variable Phi. POTEN ?s aqu? una matriu 72x36x29x30 (lonxlatxilevxtime). Compte! sembla que est? a l'inrev?s a l'especificaci? de l'arxiu netCDF quan es fa un "ncdump -c s_sep05_fism_099.nc" al linux.
lon = pi/180*netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lon'))';
lon = [lon, lon(1) + 2*pi]; % Afegeixo la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.
colat = pi/180*(90 - wrev(netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lat'))));
colat = colat(1: end/2); % M'interessa nom?s l'hemisferi Nord.
[lonMat, colatMat] = meshgrid(lon, colat); % Genero la malla.
Phi_t_ilev = NaN(size(Phi, 1), size(Phi, 2));
Phi_t_ilev(:) = Phi(:, :, Posilev, PosUT); % M'interessa el potencial per a una al?ada determinada (donada per ilev) i per a un temps UT determinat.
Phi_t_ilev = flipud(Phi_t_ilev'); % Inverteixo amb flipud, ja que la latitud comen?a a l'hemisferi Sud.
Phi_t_ilev = Phi_t_ilev(1: end/2, :); % M'interessa nom?s l'hemisferi Nord.
Phi_t_ilev = [Phi_t_ilev, Phi_t_ilev(:, 1)]; % Afegeixo la Phi corresponent a la longitud positiva coincident amb la longitud inicial + 360?. De cara a dibuixar.

% Dibuixo les dades de Phi en un mapa polar:
figure
[x_0, y_0] = pol2cart(lonMat, 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?).
h = polar(-pi: pi/100: pi, 60*ones(1, 201), '-k'); delete(h)
hold on
contourf(x_0, y_0, Phi_t_ilev, 'LineStyle', 'none', 'DisplayName', 'Phi(V)')

% Pinto meridians i paral?lels geogr?fics:
polar(-pi: pi/100: pi, 30*ones(1, 201), ':k'); polar(-pi: pi/100: pi, 60*ones(1, 201), ':k');
for i = 0:30:359
    polar([i i]*pi/180, [0 90],  ':k')
end
polar((24-UT)/12*pi*[1 1], [0 90], 'k') % Marco amb una l?nia fosca el meridi? de mitjanit.
xlim([-90 90]); ylim([-90 90])
xlabel('270', 'FontSize', 8);
ylabel('0', 'FontSize', 8);
%view(90 - 15*UT, 90)
view(90, 90)