%function plotMA

%CS: Script to plot MedsASCII files: calls readMA.m
%CS: 6 profiles to a page with all QC flags and info
%CS: Claire Spillman 24/5/06

clear all
close all
clc
mapplot=0;
%CS: Choose file to process
disp('*** PLOTTING MA FILE ***');
disp(' ')
disp('Current directory: ');
disp(pwd);
disp(' ')
disp('MA files in directory: ');
ls *.MA
disp('')
inputfile = input('MA file to plot (without .MA): ','s');
inputfile=[inputfile '.MA'];
%CS: Initialise - USER INPUTS
prof_no = 0;
mindepth = 0;
%maxdepth = 1500;
maxdepth = 1100;
mintemp = -2;
maxtemp = 32;
plotposn = [0.1  0.55 0.2 0.4;...
            0.4 0.55 0.2 0.4;...
            0.7 0.55 0.2 0.4;...
            0.1 0.05 0.2 0.4;...
            0.4 0.05 0.2 0.4;...
            0.7 0.05 0.2 0.4];
mon = {'Jan','Feb','Mar','Apr','May','Jun',...
       'Jul','Aug','Sep','Oct','Nov','Dec'};       

%%%%%%%%%%%%%% CYCLE THRU PROFILES %%%%%%%%%%%%%%%%%

%CS: Read in MA file
fid=fopen(inputfile,'r');
n=0
% Cycle through each profile
while (~feof(fid))
    n=n+1
    %CS: Increment profile number
    prof_no=prof_no+1;
    profiledata=readMA(fid,prof_no);
    
    %find out which profile is temperature:
    
    pt=strmatch('TEMP',profiledata.prof_type(:,:));
    temp=profiledata.profparm(pt,:);

    if profiledata.nprof > 1
        %setup salinity here - not used yet
    end
    lati(n)=profiledata.lat;
    longi(n)=profiledata.lon;
        
    %CS: Get number of profiles (= number of plots)
    nprofs = profiledata.nprof;
    if nprofs > 1, disp('nprofs > 1');end

    %%%%%%%%%%%%%% EXTRACT DATA %%%%%%%%%%%%%%%%%
    
    %CS Extract information for plots
    year = num2str(profiledata.year);
    month = num2str(profiledata.month);
    day = num2str(profiledata.day);
    time = profiledata.time;
    lat = profiledata.lat;
    lon = profiledata.lon;
    cruiseID = char(profiledata.cruiseID);
    
    %CS: Extract data and QC values to plot
    depth = profiledata.depth;
    depthQC = char(profiledata.depresQ);
    temp = profiledata.profparm;temp=change(temp,'>',99,NaN);
    tempQC = char(profiledata.profQparm);
    ddepth = 0;
    
    %CS: Actual QC flag i.e. TE, LE, TO etc
    nhists = profiledata.nhists;
    QCflag = cellstr(profiledata.Actcode);
    %CS: Depth for each QC flag
    QCfdep = profiledata.AuxID;
    
    %CS: Work out callsign
    surfpcode = char(profiledata.surfpcode);
    surfparm = cellstr(profiledata.surfparm);
    [m,nn]=size(surfpcode)
    callsign='unknown';
    for i = 1:m
        if surfpcode(i,1:4) == 'GCLL', callsign = char(surfparm(i)); end
    end
    
    %CS: Turn time into a string (hh:mm)
    %if time < 100, time = (time*100); end
    t = char(num2str(time));
    if length(t) == 3
        tt = [t(1),':',t(2:3)];
    elseif length(t) == 2
        tt = ['0:',t(1:2)];
    elseif length(t) == 1
        tt = ['0:0',t(1)];
    else    
        tt = [t(1:2),':',t(3:4)];
    end
    
    %%%%%%%%%%%%%%%%%% PLOTTING %%%%%%%%%%%%%%%%%%%
    
    for k = 1%:nprofs

        % Increment plot number if more than one profile
        prof_no=prof_no+(k-1);
        
        %CS: Number of plot
        nplot = rem(prof_no,6);
        if nplot == 0, nplot = 6; end
    
        %CS: need to do something here if nprof > 1 add to nplot??
        if nplot == 1
            % Create new figure for 6 new plots
            figure('position',[40 40 1000 700])
            set(gcf,'DefaultAxesFontSize',10,'DefaultTextFontAngle','normal');
            set(gcf,'PaperPositionMode','auto','color','w');
            set(gcf,'Defaulttextcolor','k','InvertHardCopy','on');
            set(gcf,'PaperOrientation','landscape');
            main_axes = gca;
            set(main_axes,'position',[0.065 0.0997 0.9500 0.800]);
            ftitle = ['CallSign: ',callsign,'    CruiseID: ',cruiseID,...
                      '    Temp (^oC)    ',mon{str2num(month)},' ',num2str(year),];
        end

        %CS: Set up subplot to plot profile    
        subplot('position',plotposn(nplot,:))
        %CS: Plot data
        plot(squeeze(temp(k,:)),squeeze(depth(k,:)),'k');
        hold on

        %CS: Plot QC flags
        p1 = [mintemp maxtemp];
        qc1 = find(tempQC(k,:) == '1');
        qc2 = find(tempQC(k,:) == '2');
        %CS: Plot dividing line where QC= 1->2
        if ~isempty(qc1)
             %CS: Deepest good depth
             ddepth = depth(max(qc1)); 
             p2 = [depth(k,max(qc1)) depth(k,max(qc1))];
             if p2(1) <= maxdepth
                plot(p1,p2,':g');
                try
                    qcnext=tempQC(k,max(qc1)+1)
                catch
                    qcnext=tempQC(k,max(qc1))
                end
                   if(qcnext~='0')
                text(0.76*maxtemp,p2(1),['QC: 1->' qcnext],'color','g','fontweight','bold',...
                     'verticalalignment','baseline');
                   end
             else
                 text(0.76*maxtemp,0.05*maxdepth,['QC: ',num2str(tempQC(10))],...
                      'color','b','fontweight','bold');
             end
        end     
        %CS: Plot dividing line where QC= 2->3
        if ~isempty(qc2)
             %CS: Update deepest good depth
             ddepth = depth(max(qc2)); 
             p2 = [depth(k,max(qc2)) depth(k,max(qc2))];
             if p2(1) <= maxdepth
                plot(p1,p2,':b');
                try
                    qcnext=tempQC(k,max(qc2)+1)
                catch
                    qcnext=tempQC(k,max(qc2))
                end
                if(qcnext~='0')
                text(0.76*maxtemp,p2(1),['QC: 2->' qcnext],'color','b','fontweight','bold', ...
                     'verticalalignment','baseline');
                   end
            end
        end 
        %CS: Put QC value of tenth value if QC > 2
        if isempty(qc1) & isempty(qc2)
            text(0.76*maxtemp,0.05*maxdepth,['QC: ',num2str(tempQC(min(10,profiledata.ndep(pt))))],...
                 'color','r','fontweight','bold');
        end

        %CS: Sort out QC flags into strings for plot
        cnt = 0; depqc = []; flgqc = [];
        for m = 1:nhists
            %CS: Count up CS flags
            if char(QCflag(m,:)) == 'CS', cnt = cnt+1; end
            %CS: Don't want more than one CS flag
            if char(QCflag(m,:)) == 'CS' &  cnt > 1 
            %CS: QC flags at same depths near the surface
            elseif m > 1 & QCfdep(m) == QCfdep(m-1) & QCfdep(m) < mindepth+2
                 flgqc(end) = cellstr([char(flgqc(end)),' ',char(QCflag(m))]);
            %CS: QC flags at same depths 
            elseif m > 1 & QCfdep(m) == QCfdep(m-1)
                 flgqc(end) = cellstr([char(flgqc(end)),' ',char(QCflag(m))]);
            %CS: Flags near the surface    
            elseif QCfdep(m) < mindepth+2
                 flgqc = [flgqc; QCflag(m)];
                 depqc = [depqc; mindepth];
            %CS: All other flags
            else
                 flgqc = [flgqc; QCflag(m)];
                 depqc = [depqc; QCfdep(m)];
            end
        end 
        % Put QC flags on plot
        for m = 1:length(depqc)
            if depqc(m) < maxdepth
                text(mintemp+0.015*(maxtemp-mintemp),depqc(m),flgqc(m),...
                      'verticalalignment','baseline');
            end      
        end   

        %CS: Put on SST (temp at 5m)
        if ddepth > 0
            m5 = find(depth >= 5);
            text(0.72*maxtemp,mindepth-0.035*(maxdepth-mindepth),...
                ['SST: ',num2str(temp(min(m5))),'C'],'color','k');
        else
            text(0.72*maxtemp,mindepth-0.035*(maxdepth-mindepth),['SST: -'],'color','k');
        end    
        hold off

        %CS: Add extra info to plot
        xlabel('Temp (^oC)');
        ylabel('Depth (m)');
        ylim([mindepth maxdepth]);
        xlim([mintemp maxtemp]);
        set(gca,'ydir','reverse');

        %CS: Add in title on plot
        date = [day '/' month '/' year '  ' tt];
        txt1 = ['Prof ',num2str(prof_no),':  ',date,'  ', ...
                num2str(lat),' ',num2str(lon)];
        xlabel(txt1);

        %CS Add max good depth to plots
        if (ddepth + 0.06*maxdepth) >= 0.95*maxdepth
            text(0.66*maxtemp,0.95*maxdepth,['Max: ',num2str(ddepth),'m']);
        elseif ddepth == 0
            %CS: Don't put on plot
        else   
            %CS Need to move below QC line
            text(0.66*maxtemp,(ddepth + 0.06*maxdepth),['Max: ',num2str(ddepth),'m']);
        end

        %CS: Figure title
        if nplot == 1 
            startprof=prof_no;
 %hh=title(ftitle,'fontsize',12,...
h2=text(-100,-500,ftitle,'fontsize',12,...
    'fontweight','bold','horizontalalignment','center',...
'Units','points','position',[-50 1]);
set(h2,'Rotation',90);
%           text(1.2*maxtemp,mindepth-0.15*(maxdepth-mindepth),ftitle,'fontsize',12,...
%'fontweight','bold','horizontalalignment','center',...
%                 'Units','points','position',[.2 .9 ])
        end 
    
        %%%%%%%%%%%%%%%%%% OUTPUT PLOTS %%%%%%%%%%%%%%%%%%%

        if nplot == 6 
            
            %CS: Can't have blanks in plot filenames
            blk = find(cruiseID == ' ');
            if ~isempty(blk), cruiseID(blk) = '-'; end

            %CS: Create filename and output to .jpg
            jpg_str = [callsign,'_',cruiseID,'_',mon{str2num(month)},'',num2str(year),'_',...
                       num2str(startprof),'-',num2str(prof_no),'.jpg']
                       jp=strfind(jpg_str,'/');
                       if(~isempty(jp));jpg_str(jp)='_';end
            screen2jpg(jpg_str)
            %CS: Create filename and output to .eps
            eps_str = [callsign,'_',cruiseID,'_',mon{str2num(month)},'',num2str(year),'_',...
                       num2str(startprof),'-',num2str(prof_no),'.eps']
                       jp=strfind(eps_str,'/');
                       if(~isempty(jp));eps_str(jp)='_';end
            screen2eps(eps_str,mapplot)
        end    
    end
end

%do last plot of incomplete pages:
if(nplot<6)
            blk = find(cruiseID == ' ');
            if ~isempty(blk), cruiseID(blk) = '-'; end

            %CS: Create filename and output to .jpg
            jpg_str = [callsign,'_',cruiseID,'_',mon{str2num(month)},'',num2str(year),'_',...
                       num2str(startprof),'-',num2str(prof_no),'.jpg']
                       jp=strfind(jpg_str,'/');
                       if(~isempty(jp));jpg_str(jp)='_';end
            screen2jpg(jpg_str)
            %CS: Create filename and output to .eps
            eps_str = [callsign,'_',cruiseID,'_',mon{str2num(month)},'',num2str(year),'_',...
                       num2str(startprof),'-',num2str(prof_no),'.eps']
                       jp=strfind(eps_str,'/');
                       if(~isempty(jp));eps_str(jp)='_';end
            screen2eps(eps_str,mapplot)
end
% Transfer plots to figures directory
% !mv *jpg dir/figures

plotmap_newreformat

msg = ['*** COMPLETED PLOTTING MA FILE ',inputfile,' ***'];
disp(' ')
disp(msg);

return