% clear memory

    %clear all ; 
    close all ;
    disp('start') ;

% start log

    logfile = [date '.log'];
    count = 0;
    while exist(['.\' logfile], 'file') == 2
        count = count + 1;
        logfile = [date '-' num2str(count) '.log'];
    end %while
    diary(logfile);

% set run parameters

    % check for parameters file and copy master if missing  
    if exist('.\setrunparameters.m', 'file') ~= 2
        try
            copyfile('d:\mje\matlab\setrunparameters.m','.')
        catch err
            disp('error 834722 copying setrunparameters.m');
        end %try
    end %if
    % run file
    try
        setrunparameters;
    catch err
        msg = 'fatal error 374285 setting run parameters';
        try 
            twit(msg);
        catch err2
        end
        disp(msg);
        rethrow(err);
    end %try
    
% combine field maps

    if combine == 1
        disp('combining field maps');
        tic;
        try
            combinefieldmaps(inputfieldmaps, fieldmapfile);
        catch err
            disp('error 274823 combining field maps');
        end %try
        toc;
        if setpauses == 1
            disp('press a key to continue...');
            pause;
        end %if
    end %if

% convert field maps to gdf

    if convert == 1
        disp('converting field map to gdf');
        tic;
        try
            gptcom(['asci2gdf -v -o fields.gdf ' fieldmapfile]);
        catch err
            disp('error 834628 converting field map to gdf');
        end %try
        toc;
        if setpauses == 1
            disp('press a key to continue...');
            pause;
        end %if
    end %if
    
% run particle tracking

    if rungpt == 1
        % check for gpt input file and copy master if missing  
        if exist('input.in', 'file') ~= 2
            try
                copyfile([gptfiles 'input.in'],'.')
            catch err
                disp('error 364829 copying input.in file');
            end %try
        end %if
        % start run
        disp('running general particle tracer');
        gptcommand = 'gpt -v -o particles.gdf input.in';
        gptcommand = [gptcommand ' npart=' num2str(npart)];
        gptcommand = [gptcommand ' Eu=' num2str(Eu)];
        gptcommand = [gptcommand ' freq=' num2str(freq)];
        gptcommand = [gptcommand ' ffac=' num2str(ffac)];
        gptcommand = [gptcommand ' rfqfac=' num2str(rfqfac)];
        gptcommand = [gptcommand ' tsim=' num2str(tsim)];
        gptcommand = [gptcommand ' tstep=' num2str(tstep)];
        gptcommand = [gptcommand ' zstep=' num2str(zstep)];
        tic;
        try
            gptcom(gptcommand);
        catch err
            disp('error 927413 running gpt');
        end %try
        toc;
        if setpauses == 1
            disp('press a key to continue...');
            pause;
        end %if
    end %if
    
% calculate trajectories

    if rungdftrans == 1
        disp('calculating trajectories');
        tic;
        try
            gptcom('gdftrans -o trajectories.gdf particles.gdf time x y z Bx By Bz rxy');
        catch err
            disp('error 936274 calculating trajectories');
        end %try
        toc;
        if setpauses == 1
            disp('press a key to continue...');
            pause;
        end %if
    end %if

% read data into matlab

    % trajmat only from trajectory file
    disp('reading trajectory file');
    tic;
    try
        [timemat, posmat, trajmat] = gdf2mat('trajectories.gdf', [], [], [], 'all') ;
    catch err
        msg = 'fatal error 742643 reading trajectory data';
        try 
            twit(msg);
        catch err2
        end
        disp(msg);
        rethrow(err);
    end %try
    toc;

    % timemat and posmat overwrite empty values from above
    disp('reading particle file') ;
    tic;
    try
        [timemat, posmat] = gdf2mat('particles.gdf', [], [], [], 'all') ;
    catch err
        msg = 'fatal error 742644 reading particle data';
        try 
            twit(msg);
        catch err2
        end
        disp(msg);
        rethrow(err);
    end %try
    toc;
    

% calculate losses
    
    disp('calculating losses');
    tic;
    try
        [firstdat, lastdat, colours, rangemat, lossbins, binmids] = gptlosses(trajmat, ranges) ;
    catch err
        disp('error 342343 calculating losses');
    end %try
    toc;

% find special particles    
    
    disp('finding special particles') ;
    tic;
    try
        rightslice = endslice;
        rightmat = timemat(rightslice) ;
        righttime = rightmat.time ;
        bunchinds = find(rightmat.z > rfqend) ;
        bunchparts = rightmat.ID(bunchinds) ;
        medEind = find(rightmat.E(bunchinds) == median(rightmat.E(bunchinds))) ;
        medpart = rightmat.ID(medEind) ;
        screenmat = posmat(rightscreen);
        survpart = length(screenmat.ID) ;
        goodeinds = find(screenmat.E > minenergy * nuc) ;
        graphinds = find(screenmat.E > mingraphenergy * nuc) ;
        goodeids = screenmat.ID(goodeinds) ;
        graphids = screenmat.ID(graphinds) ;
        goode = screenmat.E(goodeinds) ;
        graphE = screenmat.E(graphinds) ;
        survparte = length(goodeinds) ;
    catch err
        disp('error 492040 finding special particles');
    end %try
    toc;
    
% calculations

    disp('calculating results:');
    tic;
    try
        transmission = survparte/npart;
        disp([' - transmission = ' num2str(transmission*100, '%1.1f') '%']);
        meanenergy = mean(goode);
        disp([' - mean energy = ' num2str(meanenergy/nuc/1e3, '%1.0f') ' keV/u']);
        rmsenergy = std(goode);
        disp([' - energy rms = ' num2str(rmsenergy/nuc/1e3, '%1.0f') ' keV/u']);
        if calcemittance == 1
            try
                disp(' - initial x emittance');
                initialxemittance = emitcalc(timemat(1).x.*1e3, timemat(1).xp.*1e3, timemat(1).E, m);
                disp(' - initial y emittance');
                initialyemittance = emitcalc(timemat(1).y.*1e3, timemat(1).yp.*1e3, timemat(1).E, m);
                disp(' - final x emittance');
                finalxemittance = emitcalc(screenmat.x(goodeinds).*1e3, screenmat.xp(goodeinds).*1e3, screenmat.E(goodeinds), m);
                disp(' - final y emittance');
                finalyemittance = emitcalc(screenmat.y(goodeinds).*1e3, screenmat.yp(goodeinds).*1e3, screenmat.E(goodeinds), m);
            catch err2
                disp('error 823294 calculating emittance');
            end %try
        end
    catch err
        disp('error 940203 calculating results');
    end %try
    toc;
    
% build losses figures
   
    if buildlosses == 1
        disp('building loss diagrams');
        tic;
        try
            flnum = figure;
            gpttraj4plot(timemat,trajmat,colours,endslice,1,'enhance');
            figsave(flnum,'rfqlosses.jpg','screen',0);
            if setpauses == 1
                disp('press a key to continue...');
                pause ;
            end
            close(flnum);
        catch err
            disp('error 649203 building loss diagrams');
            rethrow(err);
        end %try
        toc;
    end
    
% build energy profile

    if buildenergy == 1
        disp('building energy profile');
        tic;
        try
            fenum = figure;
            figure(fenum), hold off ;
            set(fenum,'color','white') ;
            set(fenum,'Position',[50 50 600 600]) ;
            figure(fenum), hold on ;
    %        height = 0.03 * npart + 3;
            height = energymax;
            axis([0 4 0 height]) ;
            hist(graphE.*(1e-6/nuc), 100);
            xlabel('E (MeV/u)'), ylabel('n'), ...
                title('\bf rfq energy profile');
            jpgplot(fenum) ;
            hold off ;
            figsave(fenum,'rfqenergyprofile.jpg','screen','enhance');
            if setpauses == 1
                disp('press a key to continue...');
                pause ;
            end
            close(fenum);
        catch err
            disp('error 023856 building energy profile');
        end %try
        toc;
    end

% build movie files

    if buildmov == 1
    
        disp('starting movie builds') ;
        fmnum = figure;
        loopcount = [0 endslice] ;

        for j = 1:( length(loopcount) - 1 )
   
% longitudinal profile
            
            if makelongitudinal == 1
                
                disp('starting longitudinal movie') ;
                tic;
                try
                    if makemov == 1
                        mov = avifile(['rfqlongitudinal' num2str(j) '.avi'],'fps',fps,'compression','None') ;
                    end

    %                for i = (loopcount(j) + 1):loopcount(j+1)
                    for i = 1:200

                        figure(fmnum), hold off ;
                        set(fmnum,'color','white') ;
                        set(fmnum,'Position',[50 50 1000 150]) ;
                        %            medpartl = find(timemat(i).ID == medpart) ;
                        %            blen = (timemat(i).Bz(medpartl).*c)./freq ;
                        %            piclen = blen.*5 ; halflen = piclen./2 ;
                        zmin = 0; %timemat(i).z(medpartl) - halflen ;
                        zmax = rfqend + 0.05; %zmin + piclen ;

                        looppartinds = find( (timemat(i).z >= zmin) & (timemat(i).z <= zmax) ) ;

                        figure(fmnum), hold on ;
                        plot(timemat(i).z(looppartinds),timemat(i).y(looppartinds).*1e3,'xk') ;
                        xlabel('z (m)'), ylabel('y (mm)'), ...
                            title('\bf rfq longitudinal beam profile') ;
                        axis([zmin zmax -5 5]) ;
                        jpgplot(fmnum) ;
                        hold off ;

                        if makemov == 1
                            F = getframe(fmnum) ;
                            mov = addframe(mov,F) ;
                        elseif setpauses
                            disp('press a key to continue...');
                            pause ;
                        end

                        close(fmnum) ;

                    end

                    if makemov == 1
                        mov = close(mov) ;
                    end                
                catch err
                    disp('error 723953 building longitudinal movie');
                end %try                
                toc;
            
            end

% xy plane
            
            if maketransverse == 1
                
                disp('starting transverse movie');
                tic;
                try
                    if makemov == 1
                        mov = avifile(['rfqxy' num2str(j) '.avi'],'fps',fps,'compression','None') ;
                    end

    %                for i = (loopcount(j) + 1):loopcount(j+1)
                    for i = 1:200

                        figure(fmnum), hold off ;
                        set(fmnum,'color','white') ;
                        set(fmnum,'Position',[50 50 900 600]) ;
                        %            medpartl = find(timemat(i).ID == medpart) ;
                        %            blen = (timemat(i).Bz(medpartl).*c)./freq ;
                        %            piclen = blen.*5 ; halflen = piclen./2 ;
                        xmin = -3.6; %zmin = timemat(i).z(medpartl) - halflen ;
                        xmax = 3.6; %zmax = zmin + piclen ;
                        %looppartinds = find( (timemat(i).z >= zmin) & (timemat(i).z <= zmax) ) ;
                        looppartinds = find( (timemat(i).x >= xmin) & (timemat(i).x <= xmax) ) ;

                        figure(fmnum), hold on ;
                        plot(timemat(i).x(looppartinds).*1e3,timemat(i).y(looppartinds).*1e3,'xk') ;
                        xlabel('x (mm)'), ylabel('y (mm)'), ...
                            title('\bf rfq transverse beam profile') ;
                        %axis tight ;
                        %custax = axis ;
                        %custax(1) = zmin ; custax(2) = zmax ;
                        axis([-3.6 3.6 -3.6 3.6]); %axis(custax) ;
                        jpgplot(fmnum) ;
                        hold off ;

                        if makemov == 1
                            frame(i) = getframe(fmnum) ;
                            mov = addframe(mov,frame(i)) ;
                        elseif setpauses
                            disp('press a key to continue...');
                            pause ;
                        end

                        close(fmnum) ;

                    end

                    if makemov == 1
                        mov = close(mov) ;
                    end
                catch err
                    disp('error 823619 building transverse movie');
                end %try
                toc;
                
            end
            
% phase space  

            if makephasex == 1

                disp('starting x phase movie');
                tic;
                try
                    if makemov == 1
                        mov = avifile(['rfqphasex' num2str(j) '.avi'],'fps',fps,'compression','None') ;
                    end

    %               for i = (loopcount(j) + 1):loopcount(j+1)
                    for i = 1:200

                       figure(fmnum), hold off ;
                       set(fmnum,'color','white') ;
                       set(fmnum,'Position',[50 50 900 600]) ;
                       %            medpartl = find(timemat(i).ID == medpart) ;
                       %            blen = (timemat(i).Bz(medpartl).*c)./freq ;
                       %            piclen = blen.*5 ; halflen = piclen./2 ;
                       xmin = -3; %zmin = timemat(i).z(medpartl) - halflen ;
                       xmax = 3; %zmax = zmin + piclen ;
                       %looppartinds = find( (timemat(i).z >= zmin) & (timemat(i).z <= zmax) ) ;
                       looppartinds = find( (timemat(i).x >= xmin) & (timemat(i).x <= xmax) ) ;

                       figure(fmnum), hold on ;
                       plot(timemat(i).x(looppartinds).*1e3,timemat(i).xp(looppartinds).*1e3,'xk') ;
                       xlabel('x (mm)'), ylabel('x'' (mrad)'), ...
                           title('\bf rfq transverse emittance') ;
                       %axis tight ;
                       %custax = axis ;
                       %custax(1) = zmin ; custax(2) = zmax ;
                       axis([-5 5 -100 100]); %axis(custax) ;
                       jpgplot(fmnum) ;
                       hold off ;

                       if makemov == 1
                           frame(i) = getframe(fmnum) ;
                           mov = addframe(mov,frame(i)) ;
                       elseif setpauses
                           pause ;
                       end

                       close(fmnum) ;

                   end

                   if makemov == 1
                       mov = close(mov) ;
                   end
               catch err
                    disp('error 834730 building x phase movie');
               end %try
               toc;

            end    

            if makephasey == 1

                disp('starting y phase movie');
                tic;
                try
                    if makemov == 1
                        mov = avifile(['rfqphasey' num2str(j) '.avi'],'fps',fps,'compression','None') ;
                    end

    %               for i = (loopcount(j) + 1):loopcount(j+1)
                    for i = 1:200

                       figure(fmnum), hold off ;
                       set(fmnum,'color','white') ;
                       set(fmnum,'Position',[50 50 900 600]) ;
                       %            medpartl = find(timemat(i).ID == medpart) ;
                       %            blen = (timemat(i).Bz(medpartl).*c)./freq ;
                       %            piclen = blen.*5 ; halflen = piclen./2 ;
                       ymin = -3; %zmin = timemat(i).z(medpartl) - halflen ;
                       ymax = 3; %zmax = zmin + piclen ;
                       %looppartinds = find( (timemat(i).z >= zmin) & (timemat(i).z <= zmax) ) ;
                       looppartinds = find( (timemat(i).y >= ymin) & (timemat(i).y <= ymax) ) ;

                       figure(fmnum), hold on ;
                       plot(timemat(i).y(looppartinds).*1e3,timemat(i).yp(looppartinds).*1e3,'xk') ;
                       xlabel('y (mm)'), ylabel('y'' (mrad)'), ...
                           title('\bf rfq transverse emittance') ;
                       %axis tight ;
                       %custax = axis ;
                       %custax(1) = zmin ; custax(2) = zmax ;
                       axis([-5 5 -100 100]); %axis(custax) ;
                       jpgplot(fmnum) ;
                       hold off ;

                       if makemov == 1
                           frame(i) = getframe(fmnum) ;
                           mov = addframe(mov,frame(i)) ;
                       elseif setpauses
                           pause ;
                       end

                       close(fmnum) ;

                   end

                   if makemov == 1
                       mov = close(mov) ;
                   end
               
                catch err
                    disp('error 723950 building y phase movie');
                end %try
               toc;

            end
            
        end

    end

    % notify completion
    msg = ['run complete in folder ' pwd ' for user ' getenv('USERNAME')];
    try 
        twit(msg);
    catch err2
    end
    disp('end.') ;
    
    % stop log
    
    diary off;

    return

