%compares measured dispersion patterns with MAD 
clear all;
close all;


dir = '/home/skowron/dispdata'; %directory with data files
%dir = '/mnt/psop/Psop/Archives/data/CTF3/2007/week25/ResponseMtxMeas/';

homedir = '/mnt/ctfop/Public/dispmeasur'; %dir were we start and we want to end
modelhome = '/home/skowron/cern/ctfmod/cr'; %dir with the madx model

measbpms=['0130';
          '0155';
          '0195';
          '0208';
          '0248';
          '0275';
          '0305';
          '0395';
          '0425';
          '0475';
          '0570';
          '0650';
          '0695';
          '0708';
          '0775';
          '0895';
          '0925';
          '1025';
          '1070'];

[nmons, tmpi]  = size(measbpms);

fprintf(1,'Available BPMS %d\n',nmons);



%%%%% Weaker optics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   datasets= { 'BPX CR DISP  ref_1st turn orbit 1% magnets up scaled smooth optics.xls';
               'BPX CR DISP  ref_1st turn orbit 0.5% magnets up scaled smooth optics.xls';
               'BPX CR DISP  ref_1st turn orbit 0.5% magnets down scaled smooth optics.xls' };
   scalingfactors = [0.01, 0.005, -0.005];
   currentsfile='/home/skowron/ctf3emitdata/currents7332367735.dat';
%%%%% ------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Normal optics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   datasets= { 'BPX CR RC dispersion 1nd turn MKS13 20degvariation.xls';
               'BPX CR RC dispersion 2nd turn MKS13 20degvariation.xls'};
   scalingfactors = [1.0, -1.0];
   currentsfile='/home/skowron/ctf3emitdata/currents7332365434.dat';
%%%%% ------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%                                      %%%%%%%%%
%%%%%%%% CALCULATE MODEL PREVIEWED DISPERSIO  %%%%%%%%%
%%%%%%%%                                      %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%go to online model directory (containing the MADX script)
cd (modelhome);

cpycmd=['cp ' currentsfile ' currents '];
system(cpycmd);

fpo=fopen('currents','a');
fprintf(fpo,'EN=103;\n'); %energy was not stored so we need to put some
fclose(fpo);

fprintf(1,'Running MADX....\n'); %energy was not stored so we need to put some
system('madxp < dispersion.madx > out.txt'); %run the model

twissdata = importdata('dispersion.twiss',' ',47); %import the calculated orbit

cd(homedir);


ltext = length(twissdata.textdata);
[nrow,ncol] = size(twissdata.data);

fprintf(1,'Imported %d text rows, %d x %d numbers\n', ltext, nrow, ncol);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%                                      %%%%%%%%%
%%%%%%%% FIND BPM POSITIONS FROM TWISS TABLE  %%%%%%%%%
%%%%%%%%                                      %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bpmspos = zeros(nmons,1);

for m=1:nmons,
  for i=1:nrow,
     amonitorname = char(twissdata.textdata(ltext - nrow + i));   
%     fprintf(1,'Comparing %s %s\n',amonitorname,measbpms(m,:));
     if ( strfind(amonitorname,measbpms(m,:)) )
         bpmspos(m) = twissdata.data(i,1);
         fprintf(1,'Measured monitor no %d=%s is named %s and it is at pos %f\n',m,measbpms(m,:),amonitorname,bpmspos(m));
         break;
     end
  end
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%                                  %%%%%%%%%
%%%%%%%% GET THE DATA FROM EXCEL SHEETS  %%%%%%%%%
%%%%%%%%                                  %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nds,lds] = size(datasets);

fprintf(1,'datasets %d %d\n',nds, lds);

filename=' ';

disppatterns = zeros(1,1);

for ds=1:nds, %loop over available data sets

   filename=[dir '/' datasets{ds}];

   fprintf(1,'%s\n',filename);


   num = xlsread(filename, 'Diff');
   
   for i=1:length(num(:,6)),
     disppatterns(i,ds) = num(i,6)/scalingfactors(ds);
   end  

end


avgdisp = mean(disppatterns,2);

scaleatmon = 4;
scalefactor = twissdata.data(scaleatmon,2)/avgdisp(scaleatmon);

fprintf(1,'scalefactor is %f = %f / %f \n',scalefactor, twissdata.data(scaleatmon,2), avgdisp(scaleatmon));

%scalefactor = -0.0015;

disppatterns = disppatterns*scalefactor;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%                                  %%%%%%%%%
%%%%%%%% PLOT THE DATA FROM EXCEL SHEETS  %%%%%%%%%
%%%%%%%%                                  %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



for ds=1:nds, %loop over available data sets

   figure(1);
   
   h = plot(bpmspos,disppatterns(:,ds),'--o','LineWidth',2.);
   hold all;
   legend(h,'Horiz Dispersion Pattern');
   hold all;

%   figure(2);
%   h = plot(num(:,10),'-o','LineWidth',2.);
%   hold all;
%   legend(h,'Vert Dispersion Pattern')
%   hold all;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%                                      %%%%%%%%%
%%%%%%%%   PLOT MODEL PREVIEWED DISPERSION    %%%%%%%%%
%%%%%%%%                                      %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



figure(1);  
h = plot(twissdata.data(:,1),twissdata.data(:,2),'-o','LineWidth',3.);
hold all;
legend(h,'Model Dispersion');
hold all;


