% check GPS TEC binned at particular longitudes
% the script plots median TEC values; if need to plot mean TEC, change
% variable to plot to TECARRMEAN (also, change the title and plotname). 
% %%% - Larisa Goncharenko, 10/17/2011

clear
%%%%%%%%%%%%%%%%%%%%%%%%% edit here  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dirname='/hao/aim3/emery/tiegcm_dres/data/mit_climatology/'; % directory with TEC data
loncuts=[25,90,140,175,200,250,285,345]; % specify longitudes to plot
startday=datenum(2007,12,16,0,0,0); % specify first day to plot
ndays=1;  % number of days to plot
%%%%%%%%%%%%%%%%%%%%%%%%% end of edit here %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
clf
%set(gcf,'Position',[5 5 700 800],'PaperPositionMode','auto')
%set(gcf,'Position',[150 150 700 800],'PaperPositionMode','auto')
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
%
for k=1:ndays %loop over number of days
    currentday=startday+k-1;
    [year1,month1,day1,hrs1,min1,sec1]=datevec(currentday);
    clf
    for j=1:8
        clear year month day hrs mins lat lonmin lonmax tecmean tecmed dtec stderr npts
        loncut=loncuts(j);
        filename=[dirname,'long_',num2str(loncut),...
            '_',num2str(year1),'_',num2str(month1),'_',num2str(day1),'.dat'];
        if exist(filename,'file') == 2 % check if file exists
            %alldata=importdata(infile);
            [year,month,day,hrs,mins,lat,lonmin,lonmax,tecmean,tecmed,dtec,stderr,npts]=textread(filename,'%d %d %d %f %f %f %f %f %f %f %f %f %d');
            %
            dates=datenum(year,month,day,hrs,mins,0);
            datesun=unique(dates);
            datestart=datenum(year(1),month(1),day(1),hrs(1),mins(1),0);
            datesun1=datesun-datestart;
            latbin=-90:5:90;
            tecarrmean=NaN(length(datesun),length(latbin));
            tecarrmed=NaN(length(datesun),length(latbin));
            %
            for i=1:length(tecmean)
% Emery:  Toss if dtec>4 or stderr>5 TECU - not much difference for 90E 22-Nov-2007 or 16-Dec-2007
	      if dtec(i)<4 & stderr(i)<5
                idates=find(datesun == dates(i));
                ilat=find(latbin == lat(i));
                tecarrmean(idates,ilat)=tecmean(i);
                tecarrmed(idates,ilat)=tecmed(i);
              end
            end
            %
            subplot(4,2,j)
            pcolor(datesun1,latbin,tecarrmed')
            %pcolor(datesun1,latbin,tecarrmean')
            caxis([0 30]);
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
            title(['GPS TEC med, ',datestr(datestart,'dd-mmm-yyyy'),', LON=',num2str(loncut),'E']);
            %title(['GPS TEC mean, ',datestr(datestart,'dd-mmm-yyyy'),', LON=',num2str(loncut),'E']);
            set(gca,'TickDir','out');
            datetick('x','HH:MM');
            
        end
    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);
    %
    plotname=['EmeryGPSTECmed_',datestr(datestart,'dd-mmm-yyyy')];
    %plotname=['GPSTECmean_',datestr(datestart,'dd-mmm-yyyy')];
    print('-dpng',plotname);
end
