%%
%% This program is writted for comparing the TIEGCM and IRI F2 peak data 
%% Writted by J.H. lei, DEC/2005;
clear all; close all

figure;
set(gcf,'position',[ 200    38   848   723])


fpath='F:\data\storm3\storm3_weimer19\aguwb\';
%%
%% solor cycle minimum;
for iday=1:2
    %% read tgcm output;
    Fname=[fpath 'sagu-concrit2006-' num2str(346+iday) '.nc'];
    nc=netcdf(Fname,'nowrite');
    %variable=var(nc); for i=1:length(variable); disp(name(variable{i}));end

    
    time=nc{'mtime'}(:);        % (day, ut_hour, ut_min)
    UT=(time(:,1)-346-iday)*24+time(:,2)+time(:,3)/60.0;
    clear time
    % Geogetic latitude
    Lat=nc{'lat'}(:);
    % Geogetic longtitude
    lon=nc{'lon'}(:);
    height=nc{'Z'}(:);
    tgcm=nc{'NE'}(:); 
    close(nc)
    
    clear nc
    
    [UT, index]=unique(UT);
    
    height=height(index,:,:,:);
    tgcm=tgcm(index,:,:,:);
    if iday==1,   dFoF2=ones(length(Lat),length(lon),5, 2)*nan;  dTEC=dFoF2; dhmF2=dFoF2; end

    
    for i=1:4
        j=(i-1)*2+16;
        
        FoF2=ones(length(Lat),length(lon))*nan; hmF2=FoF2; TEC=FoF2;
        UTh=j;
        ip=find(UT>=UTh);
        j=ip(1);
        %% search peak density and height;
        for ii=1:length(Lat)
            for jj=1:length(lon)
                aa=fliplr(tgcm(j,:,ii,jj)); bb=fliplr(height(j,:,ii,jj));       
                cc=diff(aa); 
                for kk=1:length(cc)-1,
                    if cc(kk>=0&cc(kk+1)<=0), Kp=kk+1; break; end
                end
                X=bb(Kp-1:Kp+1);
%                Y=sqrt(aa(Kp-1:Kp+1)/1.24e4);
                Y=aa(Kp-1:Kp+1);
                [XOUX,YOUX]=HINTPO(X,Y);
                FoF2(ii,jj)=YOUX/1e6;
                TEC(ii,jj)=-trapz(bb,aa)/1e12;
                dFoF2(ii,jj,i,iday)=YOUX/1e6;
                dTEC(ii,jj,i,iday)=TEC(ii,jj);
                
                hmF2(ii,jj)=XOUX/1e5;
                dhmF2(ii,jj,i,iday)=hmF2(ii,jj);
%                 plot(aa'/1e4,bb,'r-o'); title( [num2str(j) 'UT  ' ' Lon= '  num2str(lon(jj)) '    '  num2str(FOF2(i,ii,jj)) '  ' num2str(HmF2(i,ii,jj))]); 
%                hold on;plot(1.24*YOUX^2,XOUX,'bd')
           end
        
       end

       
    end
    
end

load coast; 

%%%%%%%
for i=1:4
    j=(i-1)*2+16;

subplot(4,3,3*i-2); 
xx=[long-360 long long+360]; yy=[lat lat lat];
 plot(xx,yy,'color',[1 1 1]*0.5,'linewidth',1.5);hold on;
out=[dTEC(:,:,i,2) dTEC(:,1,i,2)]- [dTEC(:,:,i,1) dTEC(:,1,i,1)];
x=[lon; lon(1)+360]; 
xx=[x-360; x; x+360]; zz=[out out out];
pcolor(xx,Lat,zz);  caxis([-1 1]*20);  shading interp;  %colorbar 
% set(gca,'xlim',[-180 180],'xtick',[-180:60:180],'ylim',[-87 87], 'ytick',[-90:45:90],'Fontsize',11)
dg=12*15-j*15;
set(gca,'xlim',[-180+dg-0.01 180+dg],'ylim',[-80, 80],'xtick',[-180+dg:60:180+dg],'xticklabel',[],'ytick',[-80:40:90],'FontWeight','light','layer','top','TickDir','out','box','on','fontsize',11)

xx=get(gca,'xlim');
% hp=text(xx(2)*2.6/4,60,[num2str(j) ' UT'],'color','w','Fontweight','bold');
% if i==1, title(['CMIT TEC (Dec 14-Dec 13)'],'color','b','Fontsize',12,'Fontweight','bold'); end
 ylabel('Latitude (deg)','Fontsize',11);
% if i==5, xlabel('Longitude','Fontsize',12); end
if i==4,
% set(gca,'xticklabel',[-180:60:180]);
% xlabel('Degree from solar noon')
set(gca,'xticklabel',[0:4:24],'Fontsize',11);
xlabel('Local Time (hour)','Fontsize',11)
end
%if i==1,hp=text(-105,130,'Case 1'); set(hp,'fontsize',12,'color','b'); end

longm=-180+dg:5:180+dg;
load equ.mat
longin=mod(longm,360);
latm=interp1(long0,lat0,longin);
hold on; plot(longm,latm,'k:','linewidth',1.5)

 pos=get(gca,'position'); pos(2)=pos(2)+0.01*(i-1); set(gca,'position',pos);
% pos=get(gca,'position'); set(gca,'position',[pos(1)+0.15 pos(2)+0.03 pos(3:4)])

hp=title([num2str(j) ':00 UT'],'color','k','Fontsize',11);
pos=get(hp,'position');
set(hp,'position',pos+[0 -8 0])
box on

pos=get(gca,'position'); pos(3)=pos(3)*1.08;
set(gca,'position',pos)

% h=plot([-75 -75 -75],[50 35 20],'bo','markersize',4,'markerFaceColor','b','markerEdgeColor','b');


end

pos=get(gca,'position');

low=-20.; high=20.; dt=10;
pos=[ pos(1)    pos(2)-0.1 pos(3) 0.017];
axes('position',pos);
set(gca,'xlim',[1 65],'ylim',[1 65],'xtick',[],'ytick',[])
caxis([low high]);
h=colorbar('horiz') ; set(h,'position',pos)
n=(high-low)/dt;
set(h,'position',pos,'xlim',[low high],'xtick',[low:dt:high],'xticklabel',[low:dt:high],'fontsize',11)



set(gcf,'color',[1 1 1],'inverthardcopy','off');
set(gcf,'PaperPositionMode','auto');
Fname='concirt';
saveas(gcf,Fname,'tif')


