%%
%% This program is written for plotting the median or mean CCMC or TIEGCM dres TEC,NmF2,hmF2 data at 8 lonsra
%% BA Emery, December 2011 from /aim3/emery/cetic_climatology/tec8lonclim.m
clear all; close all

% mod=1-6 for IRI, CTIPE, TIE-Kp, TIE-WT, SAMI3, TIME-AMIE (ASTRA)
mod = 4;		% BAD for if12=0 from calc_med_out.m!
mod = 1;
mod = 3;

% if5deg = 1,0 for 5 deg or 2.5 deg dlat (only possible for IGS and TIE-WT)
if5deg = 1;
if (if5deg==1)
 nl36=36;
else
 nl36=71;
end

iftec = 0;
% ifmit=1,0 for mit or IGS TEC
ifmit = 1;
nt = 24;
if (iftec==1 & ifmit==0)
 nt = 12;
end
% if12 = 0 for 61-day climatology, =1,2 for 25 day climatology (hi,lo)
if12 = 0;
ifhmf2 = 1;
if iftec==1
 vartit = ' TEC '
else
 if ifhmf2==0
  vartit = 'NmF2 '
 else
  vartit = 'hmF2 '
 end
end

if mod==1
 modtit = 'IRI'
end
if mod==2
 modtit = 'CTIPe'
end
if mod==3
 modtit = 'TIE-Kp'
end
if mod==4
 modtit = 'TIE-WT'
end
if mod==5
 modtit = 'SAMI3'
end
if mod==6
 modtit = 'TIMAMI'
end
if if12==0
 modtitda = sprintf('%s 61-d ',modtit);
else
 if if12==1
 modtitda = sprintf('%s 25-d HSS ',modtit);
 else
 modtitda = sprintf('%s 25-d low Kp ',modtit);
 end
end

% Data
 loncutsE=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot
dlat = 5;		% 36*5deg=180 deg from 90S-90N		
lat5deg=-90:dlat:85;
lat2p5deg = -87.5:2.5:87.5;	

%min120 = 1:2:23;
% Make 13 values, and then duplicate last as first in arrays or otherwise plot only 11 values - not 12!
min120 = 0:2:24;
%min60=0.5:1:23.5;
min60=0:1:24;
nlon = length(min60);
if (iftec==1 & ifmit==0)
 mingps=min120;
 if (if5deg==0)
  dlat = 2.5;
  latgps=lat2p5deg;
 else
  latgps = lat5deg;
 end
else
 mingps=min60;	
 latgps = lat5deg;
end

nlon = length(mingps);	
dlon = 360/nlon;	% dlon = 360/24=15deg for hourly; dlon=15/4=2.75 for 15 min
nlat = length(latgps)

% calculate global density (and tec) and compare with Emmert at 3 altitudes
re=6371.e5;   %  earth radius in cm
dtr = pi/180.;
areatot = 0;
for ii=1:nlat
% Use midpt latitude (+dlat/2.)
 da(ii)=re*cos((latgps(ii)+dlat/2.)*dtr)*(dlon*dtr)*(re*dlat*dtr);
 areatot = areatot+da(ii)*nlon;
end

if (iftec==1)
 ratmd = NaN(length(mingps),length(latgps),8);
 pcmd = NaN(length(mingps),length(latgps),8);
else
 ratmd = NaN(length(min60),length(latgps),8);
 pcmd = NaN(length(min60),length(latgps),8);
end

% Data
        numtec=NaN(length(mingps),length(latgps),8);
        tecarrmed=NaN(length(mingps),length(latgps),8);
        numcosmic=NaN(length(min60),length(latgps),8);
	cosmicarrmed=NaN(length(min60),length(latgps),8);

figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
    clf
for j=1:8

if iftec==1
if ifmit==1
gpsloc='MIT';
 if if12==0
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/mit_gps/long_%3.3d_07325_08020_tec_mit_all.dat',loncutsE(j))
 else
  if if12==1
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/mit_gps/long_%3.3d_07325_08020_tec_mit_hi.dat',loncutsE(j))
  else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/mit_gps/long_%3.3d_07325_08020_tec_mit_lo.dat',loncutsE(j))
  end
 end
% ut_beg glat glon1 glon2 tecmed tecmean tecstd wttec averrtec numdays (max 61)
DTECin = fopen(DTECnam,'r');
datarea(j) = 0;
gdat(j) = 0;
for nh=1:nt
 for nl=1:nlat
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
  numtec(nh,nl,j) = ingps(10);
  if (ingps(10) > 0)
   if (ingps(5) > 0.) 
    tecarrmed(nh,nl,j) = ingps(5);
   else
    numtec(nh,nl,j) = 0;
   end
  end
% Calculate 'global' mean TEC from 24h of 1 glon
% Weight with area of globe and find global gdat values
	if numtec(nh,nl,j)>0
         datarea(j) = datarea(j) + da(nl);
         gdat(j) = gdat(j) + da(nl)*tecarrmed(nh,nl,j);
        end
 end	% for 1:nlat
end	% for nh=1:nt
  gdat(j) = gdat(j)/datarea(j);
else
% Read IGS TEC every 2 hours (deleted every other glat to have 5 deg glat)
gpsloc = 'IGS';
 if if12==0
  if (if5deg==1)
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs.dat',loncutsE(j))
  else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_2p5deg.dat',loncutsE(j))
  end
 else
  if if12==1
   if (if5deg==1)
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_hi.dat',loncutsE(j))
   else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_hi_2p5deg.dat',loncutsE(j))
   end
  else
   if (if5deg==1)
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_lo.dat',loncutsE(j))
   else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_lo_2p5deg.dat',loncutsE(j))
   end
  end
 end
% utmid glatmid glonmid tecmed tecmean tecstd numdays (max 61)
DTECin = fopen(DTECnam,'r');
datarea(j) = 0;
gdat(j) = 0;
for nh=1:nt
 for nl=1:nlat
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
  numtec(nh,nl,j) = ingps(7);
    tecarrmed(nh,nl,j) = ingps(4);
% Calculate 'global' mean TEC from 24h of 1 glon
% Weight with area of globe and find global gdat values
	if numtec(nh,nl,j)>0
         datarea(j) = datarea(j) + da(nl);
         gdat(j) = gdat(j) + da(nl)*tecarrmed(nh,nl,j);
        end
 end	% for 1:nlat
end	% for nh=1:nt
  gdat(j) = gdat(j)/datarea(j);
end		% end of iftec==1
else		
 if ifhmf2==0
DCOSMICnam = sprintf('/hao/aim3/emery/cetic_climatology/usu_cosmic/COSMIC_Medians_07325_08020/data_07325_08020/NmF2_%3.3d_%3.3d_07325_08020.dat',loncutsE(j)-10,loncutsE(j)+15)
 else
DCOSMICnam = sprintf('/hao/aim3/emery/cetic_climatology/usu_cosmic/COSMIC_Medians_07325_08020/data_07325_08020/hmF2_%3.3d_%3.3d_07325_08020.dat',loncutsE(j)-10,loncutsE(j)+15)
 end
DCOSMICin = fopen(DCOSMICnam,'r');
datarea(j) = 0;
gdat(j) = 0;
for nh=1:24
 for nl=1:36
 linecosmic = fgetl(DCOSMICin);
 incosmic = sscanf(linecosmic,'%f');
  numcosmic(nh,nl,j) = incosmic(10);
  if (incosmic(10) > 0) 
   cosmicarrmed(nh,nl,j) = incosmic(8);
% Weight with area of globe and find global gdat values
         datarea(j) = datarea(j) + da(nl);
         gdat(j) = gdat(j) + da(nl)*cosmicarrmed(nh,nl,j);
  end
 end	% for 1:36
end	% for nh=1:24
  gdat(j) = gdat(j)/datarea(j);
end

%% Make plots
            subplot(4,2,j)
if (iftec==1)
    gdattit = sprintf(' %5.2f',gdat(j))
    tecarrmed(nt+1,:,:) = tecarrmed(1,:,:);
            pcolor(mingps,latgps,tecarrmed(:,:,j)')
            caxis([0 30]);		% for climatology
 if if12==0
% 2 MArch CANNOT GET gpsloc in if use just title, or gdattit!
  tectit = sprintf('%s GPS TEC 61-d median, LON=%s %s',gpsloc,num2str(loncutsE(j)),gdattit); 
            title([tectit]);
 else
  if if12==1
  tectit = sprintf('%s GPS TEC 25-d median HSS, LON=%s %s',gpsloc,num2str(loncutsE(j)),gdattit);
            title([tectit]);
  else
  tectit = sprintf('%s GPS TEC 25-d median low Kp, LON=%s %s',gpsloc,num2str(loncutsE(j)),gdattit);
            title([tectit]);
  end
 end
   
else
    cosmicarrmed(nt+1,:,:) = cosmicarrmed(1,:,:);
            pcolor(min60,latgps,cosmicarrmed(:,:,j)')
 if (ifhmf2==0)
  gdattit = sprintf(' %5.2e',gdat(j))
            caxis([0 1.5e+6]);
  tectit = sprintf('COSMIC NmF2 (cm-3) 61-d median, LON=%s %s',num2str(loncutsE(j)),gdattit);
            title([tectit]);
 else
  gdattit = sprintf(' %5.1f',gdat(j))
            caxis([200 400]);
  tectit = sprintf('COSMIC hmF2 (km) 61-d median, LON=%s %s',num2str(loncutsE(j)),gdattit);
            title([tectit]);
 end
end
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
%            set(gca,'TickDir','out');
%            datetick('x','HH:MM');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);
hold off;



%%  Read Model climatology files (TEC, nmf2, and hmf2)
figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf

for j=1:8

if (mod==1)
 if if12==0
%Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_IRI.dat',loncutsE(j))		% Bad UT/LT nmf2/hmf2
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_IRI_1500km.dat',loncutsE(j))
 else
  if if12==1
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_IRI_1500km_hi.dat',loncutsE(j))
  else
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_IRI_1500km_lo.dat',loncutsE(j))
  end
 end
end
if (mod==2)
 if if12==0
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_CTIPE.dat',loncutsE(j))
 else
  if if12==1
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_CTIPE_hi.dat',loncutsE(j))
  else
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_CTIPE_lo.dat',loncutsE(j))
  end
 end
end
if mod==3
 if if12==0
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_TIE-GCM.dat',loncutsE(j))
 else
  if if12==1
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_TIE-GCM_hi.dat',loncutsE(j))
  else
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_1_TIE-GCM_lo.dat',loncutsE(j))
  end
 end
end

if mod==4
 if if12==0
   if (if5deg==1)
%Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM.dat',loncutsE(j))
% Bad GCMb for 22UT with 15 TECU in NP at 22UT?
%Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCMb.dat',loncutsE(j))
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_alld.dat',loncutsE(j))
   else
    Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_all_2p5deg.dat',loncutsE(j))
   end
 else
  if if12==1
   if (if5deg==1)
%Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_hi.dat',loncutsE(j))
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_hid.dat',loncutsE(j))
   else
    Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_hi_2p5deg.dat',loncutsE(j))
   end
  else
   if (if5deg==1)
%Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_lo.dat',loncutsE(j))
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_lod.dat',loncutsE(j))
   else
    Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/long_%3.3d_07325_08020_2_TIE-GCM_lo_2p5deg.dat',loncutsE(j))
   end
  end
 end
end

if mod==5
 if if12==0
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/sami3_output/long_%3.3d_07325_08020_Ne_sami3.dat',loncutsE(j))
 else
  if if12==1
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/sami3_output/long_%3.3d_07325_08020_tec_sami3_hi.dat',loncutsE(j))
  else
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/sami3_output/long_%3.3d_07325_08020_tec_sami3_lo.dat',loncutsE(j))
  end
 end
end

if mod==6
 if if12==0
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/astra_amie_time/long_%3.3d_07325_08020_Ne_timami.dat',loncutsE(j))
 else
  if if12==1
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/astra_amie_time/long_%3.3d_07325_08020_tec_timami_hi.dat',loncutsE(j))
  else
Dmodnam = sprintf('/hao/aim3/emery/cetic_climatology/astra_amie_time/long_%3.3d_07325_08020_tec_timami_lo.dat',loncutsE(j))
  end
 end
end

% ut_beg glat glon1 glon2 tecmed tecmean tecstd nmf2med,mean,std hmf2med,mean,std numdaystec,hmf2 (max 61)
Dmodin = fopen(Dmodnam,'r');
area = 0;
gtec(j) = 0;
gnmax(j) = 0;
ghmax(j) = 0;
nh = 0;
for mh=1:24
 nh2=fix(mh/2);
 ifuse = 1;
 if (nh2*2==mh & iftec==1 & ifmit==0)
  ifuse = 0;
 end
 if ifuse==1
  nh = nh + 1;
 end
 for nl=1:nl36
 linemod = fgetl(Dmodin);
 inmod = sscanf(linemod,'%f');
 if (ifuse==1)
  tec(nh,nl,j) = inmod(5);
  nmax(nh,nl,j) = inmod(8);
  hmax(nh,nl,j) = inmod(11);
  area = area + da(nl);
  gtec(j) = gtec(j) + da(nl)*tec(nh,nl,j);
  gnmax(j) = gnmax(j) + da(nl)*nmax(nh,nl,j);
  ghmax(j) = ghmax(j) + da(nl)*hmax(nh,nl,j);
 end	% ifuse==1
 end	% for 1:nl36
end	% for nh=1:24
gtec(j) = gtec(j)/area;
gnmax(j) = gnmax(j)/area;
ghmax(j) = ghmax(j)/area;

% Calculate ratio of model/data and percent: lat=1:2:nlat; time=0.5*(t+t+1)
     rmodat = 0;
     percentmd = 0;
     abspcmd = 0;
     ndat = 0;
     ratarea = 0;
     for nl=1:nl36
      for nh=1:nt
if (iftec==1)
	numpt = numtec(nh,nl,j);
	dat = tecarrmed(nh,nl,j);
        tie = tec(nh,nl,j);
else
	numpt = numcosmic(nh,nl,j);
	dat = cosmicarrmed(nh,nl,j);
 if ifhmf2==0
        tie = nmax(nh,nl,j);
 else
        tie = hmax(nh,nl,j);
 end
end

        if (numpt>0 & tie>0.)
	 ndat = ndat + 1;
	ratmd(nh,nl,j) = tie/dat;
        pcmd(nh,nl,j) =100.*(tie/dat-1.);
        if (tie/dat < 1.)
          pcmd(nh,nl,j) = -100.*(dat/tie-1.);
        end
% Revise pcmd for hmf2 to be in km
	if (ifhmf2==1)
         pcmd(nh,nl,j) = tie-dat;
        end
% Weight with area of globe
  	 ratarea = ratarea + da(nl); 
	 rmodat = rmodat + da(nl)*tie/dat;
         percentmd = percentmd + da(nl)*pcmd(nh,nl,j);
         abspcmd = abspcmd + da(nl)*abs(pcmd(nh,nl,j));
	end	% for numpt>0
      end	% for nh=1:nt
     end	% for nl=1:nl36
%  Divide by area of globe with numpt>0 (instead of dividing by ndat for no weighting)
	rmodat = rmodat/ratarea;
	percentmd = percentmd/ratarea;
        abspcmd = abspcmd/ratarea;
rmd8(j) = rmodat;
pcmd8(j) = percentmd;
apcmd8(j) = abspcmd;
% 53-67 COSMIC npts vs ~1150-2650 TEC GPS npts 15min 5deg bin 96*36=3356 pts total
ndat
% Put NaN in array if -1 (missing hi lats for SAMI3)
 for nl=1:nl36
  for nh=1:nt
    if (tec(nh,nl,j)<0.)
	   tec(nh,nl,j)=NaN;
	   nmf2(nh,nl,j)=NaN;
	   hmf2(nh,nl,j)=NaN;
    end
  end
 end

%% Make plots
            subplot(4,2,j)
if (iftec==1)
   tec(nt+1,:,:) = tec(1,:,:);
            pcolor(mingps,latgps,tec(:,:,j)')
            caxis([0 30]);		% for climatology
            gmodtit = sprintf(' %5.2f',gtec(j))
else
 if (ifhmf2==0)
   nmax(nt+1,:,:) = nmax(1,:,:);
            pcolor(min60,latgps,nmax(:,:,j)')
            caxis([0 1.5e+6]);
            gmodtit = sprintf(' %5.2e',gnmax(j))
 else
   hmax(nt+1,:,:) = hmax(1,:,:);
            pcolor(min60,latgps,hmax(:,:,j)')
            caxis([200 400]);
            gmodtit = sprintf(' %5.1f',ghmax(j))
 end
end
% CHANGE NAME FOR MODEL!!
           title([vartit,' ',modtitda,'median, LON=',num2str(loncutsE(j)),gmodtit]);
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
%            set(gca,'TickDir','out');
%            datetick('x','HH:MM');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    end				% end of plotting 8 longitudes

    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);

hold off;


% Plot percent of model/data
figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
for j=1:8
%% Make plots
            subplot(4,2,j)
   pcmd(nt+1,:,:) = pcmd(1,:,:);
   pctit = sprintf('%5.0f,%5.0f',pcmd8(j),apcmd8(j))
if (iftec==1)
            pcolor(mingps,latgps,pcmd(:,:,j)')
else
            pcolor(min60,latgps,pcmd(:,:,j)')
end

% CHANGE NAME FOR MODEL!!
if (ifhmf2==0)
            title([vartit,modtitda,'LON=',num2str(loncutsE(j)),' %M/D=',pctit]);
else
            pctit = sprintf('%5.1f,%5.1f',pcmd8(j),apcmd8(j))
            title([vartit,modtitda,' LON=',num2str(loncutsE(j)),' M-Dkm av,abs=',pctit]);
end

            caxis([-100 100]);
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    hc_pos(2)=0.040;
    set(hc,'Position', hc_pos);
    hold off;





