%%% This script creates input files for gfw or reads the output of gfw to
%%% plot the signal distribution.

%%% This script must be run under linux and have access to the directory
%%% where gfw is installed (this directory is set by the variable dirname
%%% below).

%%% The plots produced are stored in the directory outname.
%%% A web page with these plots is automatically created in the same
%%% directory.

%%% To create the input files set the variable create=1
%%% Once this has been done go to the gfw directory and type 
%%% gmake gfw
%%% To generate plots set create=0 (but gmake gfw must already have been
%%% run on all input files). 

%%% The machine parameters mus be set at the beginning of the while loop
%%% below and for each parameter the variation



dirname='/oxford/svn/sp_upgrade/trunk/simulation/gfw/simfiles/';
outname='/Users/delerue/travail/smith-purcell/plots/';


%%%create or plot
create=1; %create=1 to generate files,; create=0 to generate plots

machine='SOLEIL';
%machine='FACET';

%%% Create html page
if (create==0)
    fid2=fopen([ outname 'gratings_' machine '.html'],'w');
    fprintf(fid2,'<HTML>\n');
    fprintf(fid2,'<Title> %s',['Gratings studies for the ' machine ' SP experiment ']);
    fprintf(fid2,'</Title>\n');
    fprintf(fid2,'<h1> %s',['Gratings studies for the ' machine ' SP experiment  ']);
    fprintf(fid2,'</h1>\n');
    fprintf(fid2,'<table border=3>\n');
end

%% Reads GFW template
fid=fopen([ dirname 'script_gfw_template.in']);
txt=fread(fid);
fclose(fid);

changing=1;
chgmax=120;
%chgmax=4;
ichg=1;
idxchg=1;

while((changing==1)&&(ichg<chgmax)&&(ichg<=idxchg))
    idxchg=1;
    if strcmp(machine,'SOLEIL')
    %%% Set beam and grating parameters here.
    %order=1;
    grating_period=1; %mm
%    grating_period=2; %mm
%    grating_period=3; %mm
    grating_width=10; %mm
%    grating_width=40; %mm
    blaze_angle=30; % deg
    gamma=220;
    norm_emittance=10;
    position=2; %mm
    fwhmx=0.5; %mm
    fwhmy=1.1; %mm
    charge=6e8; %e-
    bunch_length_FWHM_ps=2; %ps
    shape=1;
    asymmetry=1;
    phi_step_size=1;
    number_phi_step=8;
    
    idxchg=1; % The first change is no change
    var_chg='';
    
    %%% Set parameters variation here.
    grating_period_vals=[ 0.5 1 2 3 4 5 6];
    grating_width_vals=[1.01 2 3 4 5 6 7 8 9 10 11 12 15 20 25 30 40 50 60 80 ];
    blaze_angle_vals=[10 15 17.5 20 26.75 30 36 40 45 50 ];
    position_vals=[1 2 3 4 5 6 10 20];
    fwhmx_vals=[0.3 0.5 0.7 1];
    fwhmy_vals=[0.5 1 1.1 1.5 2];
    charge_vals=[3e8 5e8 8e8 1e9 2e9 3.2e9 6.4e9 3.2e10];
    bunch_length_FWHM_ps_vals=[1 2 3 4 5];

    elseif strcmp(machine,'FACET')
                %order=1;
%    grating_period=0.05; %mm
%    grating_period=0.25; %mm
%    grating_period=0.5; %mm
%    grating_period=1.0; %mm
    grating_period=1.49; %mm
    grating_width=20; %mm
    blaze_angle=30; % deg
    gamma=45000;
    norm_emittance=15; %mm.mrad
    position=1; %mm
    fwhmx=0.2; %mm
    fwhmy=0.2; %mm
    charge=2e10; %e-
%    bunch_length_FWHM_ps=0.1; %ps
    bunch_length_FWHM_ps=0.3; %ps
    shape=1;
    asymmetry=1;
    phi_step_size=1;
    number_phi_step=8;

    grating_period_vals=[ 0.01 0.05 0.1 0.2 0.5 1 1.49 2];
    grating_width_vals=[1.01 2 3 4 5 6 7 10 15 20 25 30 40 50 60 80 ];
    blaze_angle_vals=[10 15 17.5 20 26.75 30 36 40 45 50 ];
    position_vals=[0.5 1 2 3 4 5 6 10 20];
    fwhmx_vals=[0.3 0.5 0.7 1];
    fwhmy_vals=[0.5 1 1.1 1.5 2];
    charge_vals=[5e5 1e10 2e10 3e10];
    bunch_length_FWHM_ps_vals=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0];

    end % machine
    
    
    nchg=length(grating_period_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        grating_period=grating_period_vals(ichg-idxchg);
        var_chg='grating_period';
    end %if
    idxchg=idxchg+nchg;
    
    nchg=length(grating_width_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        grating_width=grating_width_vals(ichg-idxchg);
        var_chg='grating_width';
    end %if
    idxchg=idxchg+nchg;
    
    nchg=length(blaze_angle_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        blaze_angle=blaze_angle_vals(ichg-idxchg);
        var_chg='blaze_angle';
    end %if
    idxchg=idxchg+nchg;
    
    nchg=length(position_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        position=position_vals(ichg-idxchg);
        var_chg='position';
        
    end %if
    idxchg=idxchg+nchg;
    
    nchg=length(fwhmx_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        fwhmx=fwhmx_vals(ichg-idxchg);
        var_chg='FWHM X';
    end %if
    idxchg=idxchg+nchg;
    
    nchg=length(fwhmy_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        fwhmy=fwhmy_vals(ichg-idxchg);
        var_chg='FWHM Y';
    end %if
    idxchg=idxchg+nchg;
    
    nchg=length(charge_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        charge=charge_vals(ichg-idxchg);
        var_chg='Charge';
    end %if
    idxchg=idxchg+nchg;
    
    nchg=length(bunch_length_FWHM_ps_vals);
    if ((ichg>idxchg)&&(ichg<=(idxchg+nchg)))
        bunch_length_FWHM_ps=bunch_length_FWHM_ps_vals(ichg-idxchg);
        var_chg='Bunch length';
    end %if
    idxchg=idxchg+nchg;
    
    
    
    fname=['gfw_' machine '_gperiod_' num2str(grating_period) 'mm_gwidth_' num2str(grating_width) 'mm_blaze_' num2str(blaze_angle) '_pos_' num2str(position) 'mm_fwhmx_' num2str(fwhmx) 'mm_fwhmy_' num2str(fwhmy) 'mm_charge_' num2str(charge) '_blength_' num2str(bunch_length_FWHM_ps) 'ps' ]
    
    if (create==1)
        
        outtxt=strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(strrep(char(txt'),...
            '@change_order@','n'),...
            '@grating_period@',num2str(grating_period)),...
            '@grating_width@',num2str(grating_width)),...
            '@blaze_angle@',num2str(blaze_angle)),...
            '@gamma@',num2str(gamma)),...
            '@norm_emittance@',num2str(norm_emittance)),...
            '@position@',num2str(position)),...
            '@fwhmx@',num2str(fwhmx)),...
            '@fwhmy@',num2str(fwhmy)),...
            '@charge@',num2str(charge)),...
            '@bunch_length_FWHM_ps@',num2str(bunch_length_FWHM_ps)),...
            '@shape@',num2str(shape)),...
            '@asymmetry@',num2str(asymmetry)),...
            '@phi_step_size@',num2str(phi_step_size)),...
            '@number_phi_step@',num2str(number_phi_step)),...
            '@save_dE_file@','y'),...
            '@file_dE_name@',[  'gfw_dE' '.gfwout' ]),...
            '@save_energy_file@','y'),...
            '@file_Energy_name@',[  'gfw_Energy' '.gfwout']),...
            '@quit@','q');

        
        %% Writes GFW file
        fid=fopen([ dirname fname '.gfwin'],'w');
        fwrite(fid,outtxt);
        fclose(fid);
        [ dirname fname '.gfwin']
        
    else %if not create ==1
        
        fname
        fnamerep=strrep(fname,'_',' ');
        [angle_vals,energy_vals]=plot_gratings_energy_gfw_fcn(fname,dirname,outname)
        
        
        
        %%% Check if there is a gfw31 file
        gfw31exists=fileattrib([dirname 'gfw31/' fname '_dE' '.gfwout'])
        if (gfw31exists==1)
            fname31=fname;
        else
            fname31=strrep(fname,'FACET','SOLEIL');
            gfw31exists=fileattrib([dirname 'gfw31/' fname31 '_dE' '.gfwout'])
        end %if
        if (gfw31exists==1)
            [angle_vals31,energy_vals31]=plot_gratings_energy_gfw_fcn(fname31,[dirname 'gfw31/'], [outname 'gfw31/'])
            figure(331)
            clf;
            max(energy_vals)
            max(energy_vals31)
            (max(energy_vals)/max(energy_vals31)) 
            if (((max(energy_vals)/max(energy_vals31))<20)&&((max(energy_vals31)/max(energy_vals))<20))                 
                disp('plot')
            plot(angle_vals,energy_vals,'-r','LineWidth',3)
            hold on
            grid on
            plot(angle_vals31,energy_vals31,':b','LineWidth',2)
            else
                disp('log plot')
            semilogy(angle_vals,energy_vals,'-r','LineWidth',3)
            hold on
            grid on
            semilogy(angle_vals31,energy_vals31,':b','LineWidth',2)
            end %if
            hold off
            xlabel('theta');
            ylabel('Energy/bunch (J) (25mm diameter)');
            legend('gfw33/svn88','gfw31','Location','Best');
            title({ [fnamerep ] ; 'Expected energy '});
            print('-dpng',[ outname 'compare_gfw31_33_gratings_energy_gfw_' fname '.png']);

        end %there is a gfw31 file
        
        fprintf(fid2,'<tr>\n');
        fprintf(fid2,'<td>\n');
        
        if ~isempty(var_chg)
            fprintf(fid2,'Changing variable: <b>%s</b> <BR/>\n',var_chg);
        end %if
        
        fprintf(fid2,'grating_period: %s <BR/>\n',num2str(grating_period));
        fprintf(fid2,'grating_width: %s <BR/>\n',num2str(grating_width));
        fprintf(fid2,'blaze_angle: %s <BR/>\n',num2str(blaze_angle));
        fprintf(fid2,'gamma: %s <BR/>\n',num2str(gamma));
        fprintf(fid2,'norm_emittance: %s <BR/>\n',num2str(norm_emittance));
        fprintf(fid2,'position: %s <BR/>\n',num2str(position));
        fprintf(fid2,'fwhmx: %s <BR/>\n',num2str(fwhmx));
        fprintf(fid2,'fwhmy: %s <BR/>\n',num2str(fwhmy));
        fprintf(fid2,'charge: %s (%s pC) <BR/>\n',num2str(charge), num2str(charge/1.6e7) );
        fprintf(fid2,'bunch_length_FWHM_ps: %s <BR/>\n',num2str(bunch_length_FWHM_ps));
        fprintf(fid2,'shape: %s <BR/>\n',num2str(shape));
        fprintf(fid2,'asymmetry: %s <BR/>\n',num2str(asymmetry));
        fprintf(fid2,'phi_step_size: %s <BR/>\n',num2str(phi_step_size));
        fprintf(fid2,'number_phi_step: %s <BR/>\n',num2str(number_phi_step));
        %fprintf(fid2,'save_dE_file: %s <BR/>\n','y');
        %fprintf(fid2,'file_dE_name: %s <BR/>\n',[  'gfw_dE' '.gfwout' ]);
        %fprintf(fid2,'save_energy_file: %s <BR/>\n','y');
        %fprintf(fid2,'file_Energy_name: %s <BR/>\n',[  'gfw_Energy' '.gfwout']);
        
        fprintf(fid2,'</td>\n');
        
        imgname=[ 'gratings_efficiency_gfw_' fname '.png'];
        fprintf(fid2,'<td><A href="%s"><img src="%s" width=600></A></td>\n',imgname,imgname);
        
        imgname=[ 'gratings_energy_gfw_' fname '.png'];
        fprintf(fid2,'<td><A href="%s"><img src="%s" width=600></A></td>\n',imgname,imgname);
        
                if (gfw31exists==1)
        imgname=['compare_gfw31_33_gratings_energy_gfw_' fname '.png'];
        fprintf(fid2,'<td><A href="%s"><img src="%s" width=600></A></td>\n',imgname,imgname);
                end
        fprintf(fid2,'</tr>\n');
        
        
        
        
    end %if create
    
    ichg=ichg+1
    
    %system(['ssh pplxint6.physics.ox.ac.uk "cd cauchy ; setenv  LD_LIBRARY_PATH /system/sl4/NAG/cllux08dgl/lib/ ; ./gfw31m < ./' fname '.in "'])
end %while ichg

if (create==0)
    fprintf(fid2,'</table>\n');
    fprintf(fid2,'</html>\n');
    fclose(fid2);
end

idxchg

