classdef pamtrace1

    properties (SetAccess = public)
        basefilename
        filenames
        mutants
        data
        params

        npq
    end

    methods
        function trace = pamtrace1(filename) %constructor
            if nargin < 1
                error('provide filename')
            end

            trace.basefilename  = regexp(filename, '(.*)\.asc', 'tokens');
            trace.basefilename=filename;%trace.basefilename{1}{1};
            %search current directory for all files that contain the base
            %filename in ther title
            f=ls;
            g=dir;
            kk=1;
            for k=1:length(g)
                notnumber='[\D\<]';  %to eleminate possibility that, e.g. 320 will appear if base starts with 20
                disp(g(k).name)
                ismatch1=regexp(g(k).name,[ '\<' filename]);

                ismatch2=regexp(g(k).name,[ '\D' filename]);

                switch 1
                    case length(ismatch1) %wt
                        trace.filenames.wt=g(k).name;
                        trace.data.wt=importdata(g(k).name);
                        foo=importdata(g(k).name);

                        trace.params.(mutanttype{1}{1})=getparams(g(k).name);
                        trace.data.wt=normToOne(foo.data);

                    case length(ismatch2) %mutant
                        mutanttype=regexp(g(k).name,[ '(\w*)' '_' filename], 'tokens');
                        trace.filenames.(mutanttype{1}{1})=g(k).name;
                        foo=importdata(g(k).name);

                        trace.data.(mutanttype{1}{1})=normToOne(foo.data);

                        trace.params.(mutanttype{1}{1})=getparams(g(k).name);

                end

                kk=kk+1;


            end
            trace.mutants=fields(trace.filenames);

            %calculate differences
            for k=1:length(trace.mutants)-1
                for kk=k+1:length(trace.mutants)
                    d1=trace.data.(trace.mutants{k});
                    d2=trace.data.(trace.mutants{kk});
                    l(1)  =length(d1);
                    l(2) =length(d2);
                    tstart=[1 1];
                    tend  =l;
                    [~, tfm(1)]=max(d1(:,2));
                    [~, tfm(2)]=max(d2(:,2));
                    if l(1)~=l(2)

                        l2=min(l)-max(tfm);
                        tstart=tfm-20;
                        tend=tstart+l2;
                    end

                    trace.data.([trace.mutants{k} trace.mutants{kk}])(:,1)=trace.data.(trace.mutants{k})(tstart(1):tend(1),1);
                    trace.data.([trace.mutants{k} trace.mutants{kk}])(:,2)=trace.data.(trace.mutants{kk})(tstart(2):tend(2),2)-trace.data.(trace.mutants{k})(tstart(1):tend(1),2);


                end

            end
            for k=1:length(trace.mutants)
                trace.npq.(trace.mutants{k})=getNPQ(trace,(trace.mutants{k}))
            end

            for k=1:length(trace.mutants)
                for m=k+1:length(trace.mutants)
                    trace.npq.([trace.mutants{k} trace.mutants{m}])=getNPQdiff(trace, trace.mutants{k} , trace.mutants{m} )
                end
            end
        end%pamtrace1

        function p=getNPQ(pamtrace1, mutant)

            y=pamtrace1.data.(mutant)(:,2);
            z=diff(y);



            y2=pamtrace1.data.(mutant)(:,2);
            z=diff(y);

            nmaxpeaks=60;
            [~, j]=sort(z, 'descend');
            j=j+3;
            p.tFm=pamtrace1.data.(mutant)(j(1:nmaxpeaks),1);
            [p.tFm idx]=sort(pamtrace1.data.(mutant)(j(1:nmaxpeaks),1), 'ascend')
            p.Fmprime=pamtrace1.data.(mutant)(j(idx),2);
            [p.Fm p.Fmidx]=max(p.Fmprime);
            p.Fv=p.Fm-p.Fmprime;
            p.npq=p.Fv./p.Fmprime;
        end
        function p=getNPQdiff(pamtrace1, mutant1, mutant2)
            tpoints=linspace(pamtrace1.npq.(mutant1).tFm(1),pamtrace1.npq.(mutant1).tFm(end), 50);
            m1npq=interp1(pamtrace1.npq.(mutant1).tFm,pamtrace1.npq.(mutant1).npq, tpoints, 'nearest') ;
            m2npq=interp1(pamtrace1.npq.(mutant2).tFm,pamtrace1.npq.(mutant2).npq, tpoints, 'nearest') ;
            m1Fmprime=interp1(pamtrace1.npq.(mutant1).tFm,pamtrace1.npq.(mutant1).Fmprime, tpoints, 'nearest') ;
            m2Fmprime=interp1(pamtrace1.npq.(mutant2).tFm,pamtrace1.npq.(mutant2).Fmprime, tpoints, 'nearest') ;

            p.Fmprime=m1Fmprime-m2Fmprime;
            p.tFm=tpoints;
            p.npq=m1npq-m2npq;
        end

        function plotNPQ(pamtrace1, varargin)
            plotproperty(pamtrace1,'NPQ', varargin);
        end%plot



        function plotFm(pamtrace1, varargin)
            plotproperty(pamtrace1, 'Fm', varargin)
        end%plot



        function plotproperty(pamtrace1, property,varargin )

            if length(varargin{1})<1
                mutant='wt';
            else mutant=varargin{1}{1};
            end
            if length(varargin{1})<2
                fnum=figure;
            else fnum=varargin{1}{2};
            end
            if length(varargin{1})<3
                plotstyle='bo';
            else plotstyle=varargin{1}{3};

            end
            if length(varargin{1})<4
                offset=0;
            else offset=varargin{1}{4};
            end
            if length(varargin{1})<5
                norm='';
            else norm=varargin{1}{5};
            end
            switch norm
                case ''
                    normvalue=1;
                case 'FvFm'

                    normvalue=pamtrace1.params.(mutant).FvFm;

            end
            switch property
                case 'Fm'
                    xaxis=pamtrace1.npq.(mutant).tFm;
                    yaxis=pamtrace1.npq.(mutant).Fmprime;
                case 'NPQ'

                    xaxis=pamtrace1.params.(mutant).pulsetimes;
                    yaxis=pamtrace1.params.(mutant).npq;
                case 'pamtrace1'
                    xaxis=pamtrace1.data.(mutant)(:,1);
                    yaxis=pamtrace1.data.(mutant)(:,2);
            end
            xaxis=xaxis+offset;
            figure(fnum)
            set(gca, 'fontsize', 22)

            plot(xaxis,yaxis/normvalue, plotstyle);
            grid off
            box on
          %  title(pamtrace1.basefilename)
            hold on
          %  plot(xaxis, zeros(size(xaxis)), 'k')
            xlabel('time(seconds)')

        end%plot

        function plot(pamtrace1, varargin)

            plotproperty(pamtrace1, 'pamtrace1',varargin)
        end%plot



        function plotdiff(pamtrace1, mutant,fnum, plotstyle)
            if nargin<2
                mutant='wt';
            end
            if nargin<3
                fnum=figure
            end
            if nargin<4
                plotstyle='b';
            end



            figure(fnum)
            set(gca, 'fontsize', 22)

            plot(pamtrace1.data.(mutant)(2:end,1),diff(pamtrace1.data.(mutant)(:,2)), plotstyle);
            grid off
            title(pamtrace1.basefilename(1:5))
            hold on
            plot(pamtrace1.data.(mutant)(1:end,1), zeros(size(pamtrace1.data.(mutant)(:,1))), 'k');
        end%plot





        function plotparam(pamtrace1, mutant,fnum, plotstyle, thingtoplot)
            if nargin<2
                mutant='wt';
            end
            if nargin<3
                fnum=figure
            end
            if nargin<4
                plotstyle='b';
            end
            if nargin<5
                thingtoplot='npq';
            end


            figure(fnum)
            set(gca, 'fontsize', 22)
            xaxis=pamtrace1.params.(mutant).pulsetimes;
            s=pamtrace1.params.(mutant)
            plot(xaxis,s.(thingtoplot), plotstyle);
            grid off
            title(pamtrace1.basefilename)
            hold on
            plot(xaxis, zeros(size(xaxis)), 'k')
        end%plot


    end % methods


end %classdef



function normdata= normToOne(data)
normdata=data;
normpoint=data(1,2);
normdata(:,2)=data(:,2)/normpoint;

end


function params= getparams(name)
base=regexp(name, '(.*)\.asc', 'tokens')
filename=[base{1}{1} '_params.asc'];
f=fopen(filename);
p=textscan(f, '%s %f', 'Headerlines', 1, 'MultipleDelimsAsOne', 1, 'Whitespace', ' \b\t ?');
fclose(f);
k=1;
l=2;
m=2;
n=1;
nn=1;
params.pulsetimes=importdata('pulsetimes.txt')
params.npq=zeros(size(params.pulsetimes));
for j=1:length(p{1})
    if strcmp(p{1}{j}, 'NPQ')
        params.npq(k)=p{2}(j);
        k=k+1;
    end
    if strcmp(p{1}{j}, 'qP')
        params.qP(k)=p{2}(j);
        l=l+1;
    end
    if strcmp(p{1}{j}, 'Fs')
        params.Fs(m)=p{2}(j);
        m=m+1;
    end
    if strcmp(p{1}{j}, 'Fv/Fm')
        params.FvFm(n)=p{2}(j);
       n=n+1;
    end

    if strcmp(p{1}{j}, 'Fo')
        params.Fo(nn)=p{2}(j);
       nn=nn+1;
    end

end


timeoff=1000;

toffidx=interp1( params.pulsetimes, 1:length(params.pulsetimes),timeoff, 'nearest');
npq_qEoff=params.npq(toffidx);
end