%%
%% 27 February 2013, /hao/aim3/emery/ctmt/rd_97km.m
%% Read Xiaoli Zhang's 97km data from CTMT for Dec and Sep

clear all; close all

% data(25,145,73)=(LT,glon,glat) where LT=0:1:24, glon=0:2.5:360, glat=-90:2.5:90
glon = 0:2.5:360;
glat = -90:2.5:90;
lt=0:1:24;

%for nv=1:5
for nv=2:2
fac = 1;
if (nv==1)
infile=sprintf('HME_tira_97km_tidal_u_2.5x2.5_field_09.dat');
cmax = 57;
end
if (nv==2)
infile=sprintf('HME_tira_97km_tidal_v_2.5x2.5_field_09.dat');
cmax = 54;
end
if (nv==3)
infile=sprintf('HME_tira_97km_tidal_w_2.5x2.5_field_09.dat');
cmax = 18;
fac = 100;
end
if (nv==4)
infile=sprintf('HME_tira_97km_tidal_t_2.5x2.5_field_09.dat');
cmax = 35;
end
if (nv==5)
infile=sprintf('HME_tira_97km_tidal_d_2.5x2.5_field_09.dat');
cmax = 18;
end
cmin = -cmax;
IN = fopen(infile,'r');


% 52998 lines = (5*145+1)*73-1
for il=1:73
 for jl=1:145
  for k=1:4
   line6 = fgetl(IN);
   dat6 = sscanf(line6,'%f');
   for n=1:6
    m = (k-1)*6 + n;
    data(jl,il,m) = dat6(n)*fac;
   end
  end	% for k=1:4
  line1 = fgetl(IN);
  data(jl,il,25) = sscanf(line1,'%f')*fac;
  dat00(jl,il)=data(jl,il,25);
  dat06(jl,il) = data(jl,il,7);
  dat12(jl,il) = data(jl,il,13);
  dat18(jl,il) = data(jl,il,19);
 end	% for jl=1:145
% Skip blank line except at il=73
  if (il<73)
   lineb = fgetl(IN);
  end
end	% for il=1:73


figure;
clf;

for nt=1:4
nlt = (nt-1)*6 + 1;
subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(glon,glat,dat00');
end
if (nt==2)
 contourf(glon,glat,dat06');
end
if (nt==3)
 contourf(glon,glat,dat12');
end
if (nt==4)
 contourf(glon,glat,dat18');
end
if (nv==1)
 tit = sprintf('%dLT CTMT 97km Sep Un',lt(nlt));
end
if (nv==2)
 tit = sprintf('%dLT CTMT 97km Sep Vn',lt(nlt));
end
if (nv==3)
 tit = sprintf('%dLT CTMT 97km Sep Wn',lt(nlt));
end
if (nv==4)
 tit = sprintf('%dLT CTMT 97km Sep Tn',lt(nlt));
end
if (nv==5)
 tit = sprintf('%dLT CTMT 97km Sep pc density',lt(nlt));
end
 h=title([tit]);
shading flat
caxis([cmin cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('altitude (km)','fontsize',12,'fontweight','bold');
xlabel('latitude','fontsize',12,'fontweight','bold');	

end	% for nt=1:4
hold off;	

% Convert to UT, glon vs glat from data(145glon,73glat,25LT)
nu = 0;
for iut=0:6:18
nu = nu + 1;
ut = iut;
for jl=1:145
  rlt = ut + glon(jl)/15.;
  if (rlt>=24.)
   rlt = rlt - 24.;
  end
  ltp1 = fix(rlt)+1;
  f1 = (ltp1-rlt);
  f2 = 1.-f1;
 for il=1:73
  dataut(jl,il,nu) = f1*data(jl,il,ltp1) + f2*data(jl,il,ltp1+1);
 end
end
end


figure;
clf;

for nt=1:4
nut = (nt-1)*6;
subplot(2,2,nt)
hold on;
 contourf(glon,glat,dataut(:,:,nt)');
if (nv==1)
 tit = sprintf('%dUT CTMT 97km Sep Un',nut);
 cmax = 60;
end
if (nv==2)
 tit = sprintf('%dUT CTMT 97km Sep Vn',nut);
 cmax = 60;
end
if (nv==3)
 tit = sprintf('%dUT CTMT 97km Sep Wn',nut);
end
if (nv==4)
 tit = sprintf('%dUT CTMT 97km Sep Tn',nut);
end
if (nv==5)
 tit = sprintf('%dUT CTMT 97km Sep pc density',nut);
 cmax = 15;
end
cmin = -cmax;
 h=title([tit]);
shading flat
caxis([cmin cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude (km)','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nt=1:4
hold off;	



end	% for nv=1:5
