%
% Program to calculate spurious reflections on Symmetric channel-cut monos
% Following [1] Vaclav O. Kostroun, NIM 172(1980) 243-247
% Jordi Juanhuix and I Peral, March 2007
% juanhuix@cells.es, iperal@cells.es
setconstants;

% Reproduction of the figure from [1] is done when kostroun=1
kostroun=0;

% Monochromator characteristics
si111=1;    % type of crystal, kostroun overides it
si311=0;    % type of crystal, kostroun overides it
sym_mono=1; % type of crystal, kostroun overides it

if kostroun
   % to reproduce the example shown in figure [1], it is enough to use the
   % following parameters and min=fin=4
  D = 10;                     % dist 1st xtal-screen, m. NOT USED
  X1 = 18*MILI;               % For symmetric mono: gap between xtals, negat.
  Y1 = 2.0*MILI;
  hl=20.0*MILI;               % h dimension in Vaclav O. Kostroun, NIM 172(1980) 243-247 
  L=37.0*MILI;
  alpha = [1; 1; 1];          % reference reflection
  beta = [-1; 1; 0];          % direction of Bragg axis
  gamma = [-1; -1; 2];        % orthogonal to previous ones
  asym=16.0*DEG;              % asymmetric cut
  sym_mono=0;
  maxEnergy=13.0*KILO;              % Maximum energy of the polycromatic beam (in eV)
else
    % other cases different than kostroun publication
    if and(si111==1,si311==0)
        alpha = [1; 1; 1];          % reference reflection  
        beta  = [-1; 1; 0];         % direction of Bragg axis
        gamma = [-1; -1; 2];        % orthogonal to previous ones
    end;
    if and(si111==0,si311==1)
        alpha = [3; 1; 1];          % reference reflection  
        beta  = [-1; 2; 1];         % direction of Bragg axi
        gamma = [-1; -4; 7];        % orthogonal to previous ones
    end;
    D = 10;                   % dist 1st xtal-screen (meters)
    X1 = -3*MILI;             % For symmetric mono: gap between xtals, negat.
    Y1=0.0*MILI;              % corresponds to the Y1 dimension in [1], figure 1 
    hl=50.0*MILI;             % width of the crystal (very wide, worst case scenario), 
                              % corresponds to the h dimension in  [1],
    L=200.0*MILI;             % corresponds to the L dimension in [1], figure 1
    asym=0.0*DEG;             % asymmetric cut
    maxEnergy=100.0*KILO;              % Maximum energy of the polycromatic beam (in eV)
end 
    

% X ray characteristics
a = 5.4309; %Si  -  Si(111) 3.1354161
dspacing = a/sqrt(transpose(alpha)*alpha);
dim=1;      
%Two columns Energies in keV and in Agstroms
Energy=zeros(dim,1);             % Energy in KeV
lambda=zeros(dim,1);             % wl correspond to energy values in Angstroms 
Theta=zeros(dim,1);              % Theta angles
min=8.0;                           % Minimum energy (in KeV)
fin=8.0;                           % Maximum energy (in KeV)
maxrefl = 6;                    % maximum reflection index to be searched (or plotted)

% kostroun figure publication
if kostroun==1 
    min=4;
    fin=4;
end 

for i=min:fin
  Energy(i)=i*KILO;                       % Photon Energy, keV
  lambda(i)=HCEVA/Energy(i);              % Wavelength, A
  Theta(i) = asin(lambda(i)/2/dspacing);  % Bragg angle, rad
  if kostroun 
      Theta(i)=30.0*(pi/180.0);
      lambda(i)=2*dspacing*sin(Theta(i));
      Energy(i)=HCEVA/lambda(i);
  end;
end

% Minimum distance of the extra beams from the main diffracted beam
mindisty= 3*MILI;                              
mindistz= 3*MILI;

if kostroun
    mindisty= 15*MILI;                              
    mindistz= 8.0*MILI;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% searching for close reflection for a given Theta range
for i=min:fin

  RotBragg = [sin(Theta(i))  0 cos(Theta(i));
               0           1          0;
              -cos(Theta(i)) 0 sin(Theta(i))];
          
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Reference reflection position
  h=alpha(1) ; k=alpha(2) ; l=alpha(3);
  Halpha = (alpha(1)*h + alpha(2)*k + alpha(3)*l)/ ...
            sqrt(alpha(1)^2 + alpha(3)^2 + alpha(3)^2)/sqrt(h^2 + k^2 + l^2);
  Hbeta = (beta(1)*h + beta(2)*k + beta(3)*l)/ ...
                sqrt(beta(1)^2 + beta(3)^2 + beta(3)^2)/sqrt(h^2 + k^2 + l^2);
  Hgamma = (gamma(1)*h + gamma(2)*k + gamma(3)*l)/ ...
                sqrt(gamma(1)^2 + gamma(3)^2 + gamma(3)^2)/sqrt(h^2 + k^2 + l^2);
  H = RotBragg*[Halpha; Hbeta; Hgamma];

  if sym_mono==1 
    F = X1/((2*H(1)^2-1)*sin(Theta(i)) - 2*H(1)*H(3)*cos(Theta(i))); % SYM case
  else
    % ASYM case
    F = (X1*cos(asym)+(L/2.0*cos(asym)+Y1)*sin(asym))/((2.0*H(1)*H(1)-1)*sin(Theta(i)-asym)-2.0*H(1)*H(3)*cos(Theta(i)-asym));
  end;
  
  % Point source. For extended sources change xyz distribution
  % xyz: diffraction point at 1st xtal. Lab frame
  x = 0; y = 0; z = 0;
            
  sy0 = 2*H(1)*H(2)*F + y;
  sz0 = 2*H(1)*H(3)*F + z;
            
  hold on;
  norm=sisf(h,k,l)*(GetAtomicFormFactor('Si',0.5/dspacing))^2;
  %%%%%%%%%%%%%%%%%%
           
  %%%%%%%%%%%%%%%%%%
  % searching over all possible h,k,l, or what is the same, over all possible
  % energies, assuming maxrefl=15 seems a reasonable hkl limit

   for h = 0:maxrefl
   for k = -maxrefl:maxrefl
   for l = -maxrefl:maxrefl;
      
             if check(h,k,l)==1
                 Halpha = (alpha(1)*h + alpha(2)*k + alpha(3)*l)/ ...
                  sqrt(alpha(1)^2 + alpha(3)^2 + alpha(3)^2)/sqrt(h^2 + k^2 + l^2);
                 Hbeta = (beta(1)*h + beta(2)*k + beta(3)*l)/ ...
                  sqrt(beta(1)^2 + beta(3)^2 + beta(3)^2)/sqrt(h^2 + k^2 + l^2);
                 Hgamma = (gamma(1)*h + gamma(2)*k + gamma(3)*l)/ ...
                  sqrt(gamma(1)^2 + gamma(3)^2 + gamma(3)^2)/sqrt(h^2 + k^2 + l^2);

                 H = RotBragg*[Halpha; Hbeta; Hgamma];
                  
                 if sym_mono==1 
                   F = X1/((2*H(1)^2-1)*sin(Theta(i)) - 2*H(1)*H(3)*cos(Theta(i))); % SYM case
                 else
                   % ASYM case
                   F = (X1*cos(asym)+(L/2.0*cos(asym)+Y1)*sin(asym))/((2.0*H(1)*H(1)-1)*sin(Theta(i)-asym)-2.0*H(1)*H(3)*cos(Theta(i)-asym));
                 end;
    
                 % Point source. For extended sources change xyz distribution
                 % xyz: diffraction point at 1st xtal. Lab frame
                 x = 0; y = 0; z = 0;
            
                 sx = (2*H(1)^2-1)*F + x;
                 sy = 2*H(1)*H(2)*F + y;
                 sz = 2*H(1)*H(3)*F + z;
                 
                 screen(1) = -D;
                 screen(2) = sy;
                 screen(3) = sz;
                 dist=sqrt((screen(2)/MILI-sy0/MILI)^2+(screen(3)/MILI-sz0/MILI)^2);       % distance to the reference point
                 disty=abs(screen(2)/MILI-sy0/MILI) ; distz=abs(screen(3)/MILI-sz0/MILI) ; % distance (y and z) to the reference point 
                   
                 if and(disty < mindisty/MILI,distz < mindistz/MILI)  
                   % calculating the energy of the h k l reflection
                   s0=[-1,0,0]                      ; s=s0(:)-2*(s0*H).*H;
                   mods=sqrt((s(1)^2+s(2)^2+s(3)^2)); modH=sqrt(h^2+k^2+l^2)/a;
                   TThetaf = acos(s0*s)             ; Ef=HCEVA/(2.0*(1/modH)*sin(TThetaf/2.0));
                   if Ef < maxEnergy
                       weigth=(sisf(h,k,l)*(GetAtomicFormFactor('Si',0.5*modH))^2)/norm; % "weigth" given by F_hkl^4
                       weigth=weigth^2;
                       fprintf('Reflection: %3i %3i %3i (Ef=  %6.2f  keV) at %6.2f mm from 111, weigth= %6.2f | E_mono(keV)= %d \n',h,k,l,Ef/KILO,dist,weigth,Energy(i)/KILO);
                       currpoint = plot(screen(2)/MILI,screen(3)/MILI);
                       if [h;k;l] == alpha;
                           x0 = screen(2)/MILI; y0 = screen(3)/MILI;
                           set(currpoint,'Marker','o','MarkerSize',8,'MarkerFaceColor','r','MarkerEdgeColor','r'); % multicolor: [h k l]./maxrefl
                           text(screen(2)/MILI,screen(3)/MILI,sprintf(' %3i %3i %3i',h,k,l),'FontSize',7);
                       else
                           if(check_plot(h,k,l,alpha))
                               set(currpoint,'Marker','o','MarkerSize',weigth2size(weigth),'MarkerFaceColor','b','MarkerEdgeColor','b'); % multicolor: [h k l]./maxrefl
                               text(screen(2)/MILI,screen(3)/MILI,sprintf(' %3i %3i %3i',h,k,l),'FontSize',7);
                           end;
                       end; % finishing hkl==alpha condition
                   end % finishing energy max condition
                 end; % finishing distance condition
               end;   % finishing check reflection conditions
        end;     % finishing h,k,l loop
        end;
        end;
      end;       % finishing energy loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  % Plot Options
  title({['Spurious reflections from a Si(' num2str(alpha(1))...
    num2str(alpha(2)) num2str(alpha(3)) ') monochromator']; ...
    ['Gap ' num2str(-X1/MILI) ' mm   -   E = ' num2str(Energy(i)/KILO,'%.3f') ...
    ' keV   -   d = ' num2str(dspacing,'%.4f') ' A']})
  xlabel('x (mm)');
  ylabel('y (mm)');
  axis([(sy0-mindisty)/MILI (sy0+mindisty)/MILI (sz0-mindistz)/MILI (sz0+mindistz)/MILI]); 
  
