function vtwr002( be)
%memo44256  visualization of Thomas-Wigner rotation
%
% Thomas-Wigner rotations arise in special relativity. Consider an 
% object which is accelerated by three or more non-collinear pure 
% Lorentz transformations (boosts) in such a way that it finally 
% comes to rest again. In general we will find the object not only 
% shifted, but also rotated with respect to its orginal orientation.
% The following short video illustrates this Thomas-Wigner rotation 
% by boosting a Born-rigid grid object with constant proper 
% acceleration for a finite amount of time sequentially in five 
% directions. Born-rigidity implies that the switchover between 
% one boost to the next is synchronous in the comoving inertial
% frame and asynchronous in the laboratory frame. Color-coding of 
% the grid points with the corresponding boost number (top left) 
% visualizes this asynchronicity. The clock in the top right corner 
% shows coordinate time in the laboratory frame; lengths are given 
% in light-seconds (ls).
%

% (c) G. Beyerle 2016-2017

  if nargin < 1
    % default velocity
    be = 0.7;  
  end

  doDebug = false;
%  doDebug = true
  
  doEvHorizon = false;
  assert( doEvHorizon == false); % not yet ...
  
  disp( [mfilename ': velocity beta = ' num2str( be) ' ls/s']) 
  
%  Parameters:
%  - 5 boost legs
%  - acceleration g is constant over all legs
%  - all accelerations are g > 0

  % number of time step (per boost section!)
  if doDebug
    nofTauSteps = 3
    assert( ~doVideo && ~doAniGIF)
  else
    nofTauSteps = 100;
  end

  % number of points used to mark event horizon
  if doDebug
    nofHrzPts = 3;
  else
    nofHrzPts = 51;
  end
  
  assert( local_isodd( nofHrzPts))
    
  % size of grid object in units of light-seconds
  objSideLen = 1;
 
  % square root of number of points (vertices) in grid object
  if doDebug
    nofPtsObjOneSide = 5
  else
    nofPtsObjOneSide = 21;
  end
  
  marksiz = 10;
   
  % number of (sequential) boosts
  nofLegs = 5;
  assert( nofLegs == 5)
  
  % boost acceleration
  g = 1;
  
  disp( [mfilename ': acceleration g = ' num2str( g) ' ls/s^2'])
  
  % calculate
  % - normalized boost directions eNrm
  % - Lorentz transformation matrices LT
  % - Thomas-Wigner rotation angle thomasAngleRad
  [eNrm, allLT, LT, ...
    thomasAngleRad] = local_boostdirsandaccelerations( be, nofLegs, g);  
    
  disp( [mfilename ': Thomas-Wigner angle = ' ...
    num2str( thomasAngleRad*180/pi) ' deg'])
    
  % plot boost directions
  if false
    figNo = 1;  
    local_plotdirections( figNo, eNrm, be);
  end
       
  % initialize grid object  
  [gridObj, nofPtsObj] = local_initializegridobject( objSideLen, ...
      nofPtsObjOneSide);
  
  % traceObj_4v_fr1 : 4 x nofObjPts x nofTimeSteps
  % traceObj_4v_fr1(:,1,:)  denotes the reference point
  % timProper : proper time of all grid vertices
  % 4-positions of event horizons
  % traceHrz_4v_fr1  = 4 x (nofLegs*nofHrzPts) x nofLegs*nofTau
  % traceHrz2_4v_fr1 = 4 x nofHrzPts x nofLegs x nofLegs*nofTau
  % h : nof points on event horizon
  % t : nof time step (in each leg)
  % l : boost leg number (l = 1,...,5)
  % traceHrz_4v_fr1(:,(1:h),(1:t))         : 1st horizon on first boost
  % traceHrz_4v_fr1(:,(1:h),(1:t)+(l-1)*t) : 1st horizon on boost l>1
  
  % trace grid object 
  [traceObj_4v_fr1, timProper] = local_traceobject( nofLegs, ...
    nofTauSteps, be, eNrm, g, LT, gridObj);  
  
  % event horizon not yet functional
  traceHrz_4v_fr1 = [];
    
  % get info on corresponding boost leg number
  boostLegNo_fr1 = local_getboostlegnumber( nofPtsObjOneSide, ...
    nofTauSteps, nofLegs);
  
  % interpolate from proper times to to coordinate (lab) time
  [traceObj_4v_fr1_eqlTime, traceHrz_4v_fr1_eqlTime, timProper_eqlTime, ...
    boostLegNo_fr1_eqlTime, timLab_fr1_eqlTime] = ...
    local_interpolatetolabtime( doEvHorizon, nofHrzPts, ...
      nofLegs, nofPtsObj, traceObj_4v_fr1, traceHrz_4v_fr1, ...
      timProper, boostLegNo_fr1);   
    
% ---------------
%
%  plot grid object on its trajectory
%
    
%   maxCol = 16;
%   tmp    = parula(maxCol);
%   colorList = tmp(round( linspace( 1, maxCol, nofLegs)),:);
 
  colorList = [0, 0, 1; ...
               0, 1, 0; ...
               0, 1, 1; ...
               1, 0, 1; ...
               0, 0, 1];
    
  nofTimeSteps = size( traceObj_4v_fr1, 3);
    
  textStr = { ...
'Thomas-Wigner rotations arise in special relativity. Consider a ' ...
'Born-rigid object. It accelerates in three or more non-collinear ' ...
'pure Lorentz transformations (boosts) in such a way that it ' ...
'finally comes to rest at the starting point again. It turns out, ' ...
'that the object is not only shifted, but also rotated with respect '...
'to its orginal orientation. The visualization illustrates this ' ...
'Thomas-Wigner rotation by sequentially boosting a Born-rigid grid ' ...
'with constant proper acceleration for a finite amount of time in ' ...
'five directions. Born-rigidity implies that the switchover between ' ...
'one boost and the next is synchronous in the comoving, inertial ' ...
'reference frame and thus asynchronous in the laboratory frame. ' ...
'Color-coding of the grid points with the corresponding boost ' ...
'number visualizes this asynchronicity. The clock in the top right ' ...
'corner shows coordinate time in the laboratory frame; lengths are ' ...
'given in light-seconds (ls).', ...
   };

  disp( ' ')
  for k = 1:length(textStr)
    disp( textStr{k})
  end
  disp( ' ')
 
  
  if ~doEvHorizon
    
    tmp = traceObj_4v_fr1(2,:,:);
    minX = ( min( tmp(:)));
    maxX = ( max( tmp(:)));

    tmp = traceObj_4v_fr1(3,:,:);
    minY = ( min( tmp(:)));
    maxY = ( max( tmp(:)));

    minX = floor( minX*5) / 5 - 1e-8;
    maxX = ceil(  maxX*5) / 5 + 1e-8;
    minY = floor( minY*5) / 5 - 1e-8;
    maxY = ceil(  maxY*5) / 5 + 1e-8;

  else
    
    minX = -1.5
    maxX = +4
    minY = -2
    maxY = +4
    
  end
  
  % loop over all time step (all five legs)
  for t = 1:nofTimeSteps,

    figure(1)
    clf
    plot(NaN,NaN)
    hold on

    yOfs = 0;
    for l = 1:nofLegs,    
      sel = find( boostLegNo_fr1_eqlTime(:,t) == l);
      if isempty( sel)
        % skip if this leg not present/active at current time step
        continue
      end
      
      h = plot( traceObj_4v_fr1_eqlTime(2,sel,t), ...
                traceObj_4v_fr1_eqlTime(3,sel,t), '.', ...
                'color', colorList(l,:));
      set( h, 'markersize', marksiz)
      
      % plot boost number information
      text( 0.02, 0.95-yOfs, ...
        ['boost no. ' num2str(l)], ...
        'units', 'normalized', 'horizontalalign', 'left', ...
        'fontsize', 12, 'color', colorList(l,:))
      yOfs = yOfs+0.05;
      
    end

    h = plot( traceObj_4v_fr1_eqlTime(2,1,t), ...
              traceObj_4v_fr1_eqlTime(3,1,t), 'r.');
    set( h, 'markersize', marksiz)

    if doEvHorizon
    % plot event horizon
      for l = 1:nofLegs,    
        h = plot( traceHrz_4v_fr1_eqlTime(2,:,l,t), ...
                  traceHrz_4v_fr1_eqlTime(3,:,l,t), 'k.');
        set( h, 'markersize', marksiz)
      end
    end    
  
    % coordinate time in laboratory frame
    text( 0.98, 0.95, ...
      sprintf( '%.3f s', traceObj_4v_fr1_eqlTime(1,1,t)), ...
      'units', 'normalized', 'horizontalalign', 'right', ...
      'fontsize', 12)

    % coordinate grid
    plot( [minX,maxX], [0, 0], 'k:')
    plot( [0, 0], [minY, maxY], 'k:')
    set( h, 'markersize', 10)
    set( gca, 'xlim', [minX,maxX])
    set( gca, 'ylim', [minY,maxY])      
    set( gca, 'tickdir', 'out')      
    title( 'Visualization Thomas-Wigner rotation')
    ylabel( 'y [ls]')
    xlabel( 'x [ls]')
    daspect([1,1,1])
    
    drawnow
    
  end % --- for t = 1:nofTimeSteps, ---

  return  
    
end

% -------------------------------------------------------------------------
%
% calculate
% - normalized boost directions eNrm
% - Lorentz transformation matrices LT
% - Thomas-Wigner rotation angle thomasAngleRad
%
function [eNrmOut, allLT, LT, ...
  thomasAngleRad] = local_boostdirsandaccelerations( be, ...
  nofLegs, acc)

  % boost angles between boost direction eNrm1 and eNrm2,
  % eNrm2 and eNrm3, eNrm1 and eNrm3.
  [C12, C23, C13] = local_relrotcosangles5seg( be);
  alphaRadListOut = acos([ C12, C23, C23, C12]);
  
  eNrmOut = local_getunitvectors( alphaRadListOut);

  allLT = eye( 4);  
  for i = 1:nofLegs,
    allLT = local_lorentztrafo( eNrmOut(:,i), be) * allLT;
  end
  
  
  [ret, allPosAtEv_4v, LT] = local_minfun( alphaRadListOut, ...
    be, nofLegs, acc);
     
  allPosAtEv_4v_fr1 = zeros( size( allPosAtEv_4v)) * NaN;

  invLT = eye(4);
  
  for k = 1:nofLegs,
    
%    invLT = invLT * inv( LT(:,:,k));
    invLT = invLT / LT(:,:,k);
    
    allPosAtEv_4v_fr1(:,k) = invLT * allPosAtEv_4v(:,k);
    
  end
 
  % estimate Thomas-Wigner rotation angle
  [v, d] = eig( allLT(2:4,2:4));
  eigVal = diag( d);
  eigVal(abs( imag( eigVal)) < 1e-6) = [];
  
  if ~isempty( eigVal)
    thomasAngleRad = angle( eigVal(1));
  else    
    thomasAngleRad = NaN;
  end
 
end

% -------------------------------------------------------------------------
%
% calculate
% - trajectories
% - Lorentz transformation matrices
%
function [ret, allPosAtEv_4v, LT] = local_minfun( x, be, nofLegs, acc)

  alphaRadList = x;
    
  g = acc;
    
  alpha = [0, cumsum( alphaRadList)];
  
  eNrm = zeros( 3, nofLegs);  
  
  eNrm(1,:) = cos( alpha);
  eNrm(2,:) = sin( alpha);
    
  LT = zeros(4,4,nofLegs);

  allLT = eye(4);
  
  for k = 1:nofLegs,
    
    % LT : Lorentz trafo from frame (k) to frame (k+1)
    LT(:,:,k) = local_lorentztrafo( eNrm(:,k), be);
    allLT = LT(:,:,k) * allLT;
    
  end
  
  % boost distance
  w = atanh( be);
  sinhW = sinh( w);
  coshW = cosh( w);
  
  if length( g) == 1
    g = g * ones(1,nofLegs);
  end
  
  x = 1 ./ g .* ( coshW - 1);
  t = 1 ./ g * sinhW;
  
  posAtEv = zeros(4,1);
  
  allPosAtEv_4v = zeros(4,nofLegs);
  
  for k = 1:nofLegs,

    posAtEv = LT(:,:,k) * (posAtEv + [t(k); x(k)*eNrm(:,k)]);
    
    allPosAtEv_4v(:,k) = posAtEv;
    
  end
    
  ret = sum( allLT(1,2:4).^2) + ...
        sum( allLT(2:4,1).^2) + ...
        sum( posAtEv(2:4).^2);
      
%   ret = sum( allLT(1,2:4).^2) + ...
%         sum( allLT(2:4,1).^2);
  
end


% -------------------------------------------------------------------------
%
% plot normalized boost directions
%
function local_plotdirections( figNo, eNrm, be)

  N = size( eNrm, 2);

  figure(figNo)
  clf
  plot3(NaN,NaN,NaN)
  axis( gca);
  hold on
%   set( gca,'xlim', [-2, 10])
%   set( gca,'ylim', [-2, 10])
%   set( gca,'zlim', [-2, 10])
  daspect([1 1 1]);

  ofs = [0; 0; 0];
  
  for k = 1:N    
%    arrow( ofs, ofs + eNrm(:,k), 'facecolor', 'blue', 'edgecolor', 'blue')
    h = quiver3( ofs(1), ofs(2), ofs(3), eNrm(1,k), eNrm(2,k), eNrm(3,k));
    h.MarkerSize      = 10;
    h.LineWidth       = 2;
    h.Color           = 'blue';
    h.MaxHeadSize     = 0.3;
    h.AutoScaleFactor = 1;
%    h
%    keyboard
    ofs = ofs + eNrm(:,k);
  end
  
  title( ['\beta = ' num2str(be)])  
  view(0,90);    
  drawnow

  keyboard
  
end

% -------------------------------------------------------------------------
%
%  Lorentz trafo in direction eNrm
%
function LT = local_lorentztrafo( eNrm, be)

  ga = 1 / sqrt(1 - be^2);

  eNrmT = transpose( eNrm);
  
%   LT = eye(4) + ...
%        [(ga-1),       +ga*be*eNrmT; ...
%         +ga*be*eNrm, (ga-1)*(eNrm*eNrmT)];
  LT = eye(4) + ...
       [(ga-1),       -ga*be*eNrmT; ...
        -ga*be*eNrm, (ga-1)*(eNrm*eNrmT)];
      
end

% -------------------------------------------------------------------------
%
% from boost angles derive boost unit vectors
%
function eNrm = local_getunitvectors( alphaRad)
  
  nofLegs = length( alphaRad)+1;
  
%   alphaRad = 2 * pi / nofLegs;
%   alphaRadList = (0:nofLegs-1) * alphaRad;

  alphaRadList = [0, cumsum( alphaRad)];

  eNrm = zeros( 3, nofLegs);  
  eNrm(1,:) = cos(alphaRadList);
  eNrm(2,:) = sin(alphaRadList);

end

% -------------------------------------------------------------------------
%
%
%
%
function [crsDir1, crsDir2] = local_findperpnormals( boostDir)

  tmp = cross( [0;0;1], boostDir);
  if max( abs( tmp)) < 1e-3,
    tmp = cross( [0;1;0], boostDir);
  end
  
  crsDir1 = tmp / sqrt( sum( tmp.^2, 1));
  crsDir2 = cross( crsDir1, boostDir);

end

% -------------------------------------------------------------------------
%
%  get number of boost for each grid point
%
function boostLegNo_fr1 = local_getboostlegnumber( nofPtsObjOneSide, ...
    nofTauSteps, nofLegs)
  
  boostLegNo_fr1 = zeros( 1+nofPtsObjOneSide^2, ...
    nofTauSteps*nofLegs) * NaN;

  sel = (1:nofTauSteps);
  for i = 1:nofLegs
    boostLegNo_fr1(:,sel) = i;
    sel = sel + nofTauSteps;
  end
  
end

% -------------------------------------------------------------------------
%
%  interpolate to same coordinate time
%
function [traceObj_4v_fr1_eqlTime, traceHrz_4v_fr1_eqlTime, ...
    timProper_eqlTime, boostLegNo_fr1_eqlTime, ...
    timLab_fr1_eqlTime] = local_interpolatetolabtime( doEvHorizon, ...
    nofHrzPts, ...
    nofLegs, nofPtsObj, traceObj_4v_fr1, traceHrz_4v_fr1, ...
    timProper, boostLegNo_fr1)

  % interpolate all points to reference point time
  % (reference point is point no 1.)
  traceObj_4v_fr1_eqlTime = zeros( size( traceObj_4v_fr1)) * NaN;

  % event horizon
  traceHrz_4v_fr1_eqlTime = zeros( size( traceHrz_4v_fr1)) * NaN;

  % proper time of individual object points
  timProper_eqlTime       = zeros( size( timProper)) * NaN;

  boostLegNo_fr1_eqlTime  = zeros( size( boostLegNo_fr1)) * NaN;

  % copy reference point
  % size( traceObj_4v_fr1) = 4 x nofPtrGridObj x nofTimeSteps

  % copy grid object's reference point trajectory
  traceObj_4v_fr1_eqlTime(:,1,:) = traceObj_4v_fr1(:,1,:);

  % not  traceHrz_4v_fr1(:,1,:)  !!!
  for l = 1:nofLegs
    traceHrz_4v_fr1_eqlTime(1,1,l,:) = traceObj_4v_fr1(1,1,:);
  end
  
  % object four vector in lab frame (frame [1])

  % copy reference point's time to all other grid points
  traceObj_4v_fr1_eqlTime(:,2:end,:) = repmat( ...
      traceObj_4v_fr1_eqlTime(:,1,:), [1, nofPtsObj-1, 1]);
    
%   traceHrz_4v_fr1_eqlTime(:,2:end,:) = repmat( ...
%       traceHrz_4v_fr1_eqlTime(:,1,:), [1, nofHrzPts*nofLegs-1, 1]);
%   traceHrz_4v_fr1_eqlTime(:,2:end,:) = repmat( ...
%       traceHrz_4v_fr1_eqlTime(:,1,:), [1, nofHrzPts*nofLegs-1, 1]);

  % loop over all points of grid object
  % skip first point (reference point)
  for p = 2:nofPtsObj,

    % loop over x and y coordinate (motion is restricted to x-y plane)
    % linear interpolation of spatial coordinates
    for posId = (2:4),
      traceObj_4v_fr1_eqlTime(posId,p,:) = interp1q( ...
          squeeze( traceObj_4v_fr1(1,p,:)), ...
          squeeze( traceObj_4v_fr1(posId,p,:)), ...
          squeeze( traceObj_4v_fr1_eqlTime(1,1,:)));
    end

    timProper_eqlTime(p,:) = interp1q( ...
        squeeze( traceObj_4v_fr1(1,p,:)), ...
        transpose( timProper(p,:)), ...
        squeeze( traceObj_4v_fr1_eqlTime(1,1,:)));

    boostLegNo_fr1_eqlTime(p,:) = interp1q( ...
        squeeze( traceObj_4v_fr1(1,p,:)), ...
        transpose( boostLegNo_fr1(p,:)), ...
        squeeze( traceObj_4v_fr1_eqlTime(1,1,:)));

  end % --- for p = 2:nofPtsObj, ---

  boostLegNo_fr1_eqlTime = floor( boostLegNo_fr1_eqlTime);

  % at end of interpolation range there might be NaNs ...
  boostLegNo_fr1_eqlTime(:,end) = boostLegNo_fr1_eqlTime(:,end-1);

  
  if doEvHorizon
    
    % interpolate event horizons to same time grid
    for l = 1:nofLegs,
      for p = 1:nofHrzPts,
        for r = (2:4),

  %         traceHrz_4v_fr1_eqlTime(r,p,l,:) = interp1q( ...
  %             squeeze( traceHrz_4v_fr1(1,p,l,:)), ...
  %             squeeze( traceHrz_4v_fr1(r,p,l,:)), ...
  %             squeeze( traceHrz_4v_fr1_eqlTime(1,1,l,:)));

          x = squeeze( traceHrz_4v_fr1(1,p,l,:));
          y = squeeze( traceHrz_4v_fr1(r,p,l,:));
          xx = squeeze( traceHrz_4v_fr1_eqlTime(1,1,l,:));

          sel = ( ~isnan( x) & ~isnan( y));

          traceHrz_4v_fr1_eqlTime(r,p,l,:) = interp1( ...
              x(sel), y(sel), xx, 'previous');

        end
      end
    end
    
  end

  % interpolation of last time step may fail, so we just copy it
  tmp = traceObj_4v_fr1_eqlTime(1,:,end);
  delta = max( abs( tmp - tmp(1)));
  
  assert( delta < 1e-9)
  
  traceObj_4v_fr1_eqlTime(:,:,end) = traceObj_4v_fr1(:,:,end);
  timProper_eqlTime(:,end)         = timProper(:,end);

  timLab_fr1_eqlTime = squeeze( traceObj_4v_fr1_eqlTime(1,:,:));

end

% -------------------------------------------------------------------------
%
%  initialize grid object
%
function [gridObj, nofPtsObj] = local_initializegridobject( ...
  objSideLen, nofPtsObjOneSide)

  % grid object
  [x, y] = meshgrid( ...
             linspace( -objSideLen/2, objSideLen/2, nofPtsObjOneSide), ...
             linspace( -objSideLen/2, objSideLen/2, nofPtsObjOneSide));
  
  nofPtsObj = numel(x) + 1;
  gridObj = zeros( 3, nofPtsObj);
  gridObj(1,2:end) = x(:);
  gridObj(2,2:end) = y(:);
  
end

% -------------------------------------------------------------------------
%
% trace object (grid consisting of 'nofPts' individual points) 
%   and event horizon thru all five boost segments
%
%  traceObj_4v_fr1  : 4-D trace of (square) object in frame 1 (lab frame)
%  traceHrz_4v_fr1  : event horizon in frame 1 (lab frame)
%  traceHrz2_4v_fr1 : event horizon in frame 1 (lab frame)
%  timProper        : proper time
%
%  nofLegs     : number of boost segments
%  nofTauSteps : 
%  be          : max. velocity
%  eNrm        : boost direction in n'th frame
%  accRef      : boost acceleration of reference point
%                trailing points experience larger,
%                leading points experience smaller acceleration
%  LT          : Lorentz transformation matrices 
%  LT(:,:,k)   : Lorentz trafo from frame [k] to frame [k+1]
%  gridObj     : grid points of boosted object
%
function [traceObj_4v_fr1, timProper] = ...
  local_traceobject( nofLegs, nofTauSteps, be, eNrm, accRef, ...
  LT, gridObj) 

  % boost (proper) time of reference point
  tau = atanh( be) / accRef;
  
  % object proper time at reference point
  tauList = linspace( 0, tau, nofTauSteps);
  
  nofPtrGridObj = size( gridObj, 2);
  
  nofTimeSteps = nofLegs * nofTauSteps;
  
  % object four vector in lab frame (frame [1])
  traceObj_4v_fr1 = zeros( 4, nofPtrGridObj, nofTimeSteps) * NaN;
  
  % proper time at each object point
  timProper = zeros( nofPtrGridObj, nofTimeSteps) * NaN;
  
  % start at origin and time zero
  posGridObjPts_4v    = zeros( 4, nofPtrGridObj);
  newPosGridObjPts_4v = zeros( 4, nofPtrGridObj) * NaN;

  % initialize
  posGridObjPts_4v(2:4,:) = gridObj;

%  nofPtsEvHrzPos    = nofHrzPts * nofLegs;
%  posEvHrzPts_4v    = zeros( 4, nofPtsEvHrzPos) * NaN;

  
%  newPosEvHrzPts_4v = zeros( 4, nofPtsEvHrzPos) * NaN;
    
  %  posRef_4v = zeros( 4, 1);
  
  % loop ober all boost legs
  
  % lorentz transform from frame (k) to lab frame (frame (1))
  LTToFr1 = eye(4);
  
%   % event horizons
%   % length of horizon (of course should be infinite)
%   hrzLength = 2;
%   % number of points used to characterize the event horizon
%   % nofHrzPts = 20;
  
%   % for each time step (nofTimeSteps) we have  
%   % nofLegs*nofHrzPts  4-vector coordinates
%   traceHrzOld_4v_fr1 = zeros( 4, nofLegs*nofHrzPts, nofTimeSteps) * NaN;
%   traceHrz_4v_fr1    = zeros( 4, nofHrzPts, nofLegs, nofTimeSteps) * NaN;
    
  % loop over all five boost segments
  % (we do 'nofTauSteps' time steps at once)
  for k = (1:nofLegs),

%    disp([mfilename ': leg = ' num2str(k)])
    
    boostDir           = eNrm(:,k);
    [crsDir1, crsDir2] = local_findperpnormals( boostDir);
      
    assert( abs( norm( crsDir1) - 1) < 1e-9)
    assert( abs( norm( crsDir2) - 1) < 1e-9)
    
    % in each leg loop over all time steps
%    for t = 1:nofTauSteps,

    % first point  posGridObjPts_4v(:,1)  is reference point!
    % (current reference frame)
    posRef_m = repmat( posGridObjPts_4v(2:4,1), 1, nofPtrGridObj);

    % offset w.r.t. to reference point
    ofsA2X = posGridObjPts_4v(2:4,:) - posRef_m; 
    assert( all( ofsA2X(:,1) == 0))

%     Xcrs1  = transpose( crsDir1)  * OfsA2X;
%     Xcrs2  = transpose( crsDir2)  * OfsA2X;

    % 'boostCoord' : offset to ref point projected on (normal) boost vector     
    %   minimum value of 'boostCoord' is  -1 / accRef
    boostCoord  = transpose( boostDir) * ofsA2X;
    
    % boost factor
    L = 1 + accRef * boostCoord;
        
%    disp( [mfilename ': L = ' num2str( min(L(:))) ' ... ' num2str( max(L(:)))])
    
    % remove point beyont event horizon (L<0)
    idx = (L < 0);
    L(idx) = NaN;

    tSel = (1:nofTauSteps) + (k-1)*nofTauSteps;
            
    %  posGridObjPts_4v(2:4,1)  :  reference point

    % loop over all points of grid object, boost each one of them
    for p = (1:nofPtrGridObj),
            
      % tmpPos : 3-position in corresponding boost frame (frames [1] to [5])
      tmpPos = repmat( posGridObjPts_4v(2:4,p), 1, nofTauSteps) ...
               + boostDir / accRef * L(p) * (cosh( accRef * tauList) - 1);
                            
      % tmpTim : coordinate time of corresponding boost frame
      tmpTim = posGridObjPts_4v(1,p) + L(p) / accRef * sinh( accRef * tauList);
      
      % Lorentz trafo from current frame to lab frame (frame [1])
      traceObj_4v_fr1(:,p,tSel) = LTToFr1 * [tmpTim; tmpPos];
      
      % LT(:,:,k) : Lorentz trafo from frame (k) to frame (k+1) 
      newPosGridObjPts_4v(:,p) = LT(:,:,k) * [tmpTim(1,end); tmpPos(:,end)];
      
      % object proper time
      % timProper(p,tSel) = tauList + (k-1) * tau * L(p);
      if tSel(1) == 1,
        timProper(p,tSel) = tauList  * L(p);
      else
        timProper(p,tSel) = timProper(p,tSel(1)-1) + tauList  * L(p);
      end
          
    end  % --- for p = 1:nofPtrGridObj, ---
    
    % make sure coordinate time are all identical
    assert( max( abs( newPosGridObjPts_4v(1,:) ...
      - newPosGridObjPts_4v(1,1))) < 1e-10)
    
       
    % newPosGridObjPts_4v : 4-positions of grid points in next frame
    posGridObjPts_4v = newPosGridObjPts_4v;
    
%    posEvHrzPts_4v   = newPosEvHrzPts_4v;
    
    % LTToFr1 = LTToFr1 * inv( LT(:,:,k));
    % this numerically better than  
    %   LTToFr1 := LTToFr1 * inv( LT(:,:,k))
    LTToFr1 = LTToFr1 / LT(:,:,k);
    
  end % --- for k = 1:nofLegs, ---

  return
        
end


% -------------------------------------------------------------------------
%
% trace event horizons thru all five boost segments
%
%  traceObj_4v_fr1  : 4-D trace of (square) object in frame 1 (lab frame)
%  traceHrz_4v_fr1  : event horizon in frame 1 (lab frame)
%  timProper        : proper time
%
%  nofLegs     : number of boost segments
%  nofTauSteps : 
%  be          : max. velocity
%  eNrm        : boost direction in n'th frame
%  accRef      : boost acceleration of reference point
%                trailing points experience larger,
%                leading points experience smaller acceleration
%  LT          : Lorentz transformation matrices 
%  LT(:,:,k)   : Lorentz trafo from frame [k] to frame [k+1]
%
  
function traceHrz_4v_fr1 = local_tracehorizons( nofLegs, nofTauSteps, ...
  nofHrzPts, be, eNrm, accRef, LT, traceObj_4v_fr1, timProper, hrzLength) 

  % boost time of reference point
  tau = atanh( be) / accRef;
  
%   % object proper time at reference point
%   tauList = linspace( 0, tau, nofTauSteps);
  
%  nofPtrGridObj = size( gridObj, 2);
  
  nofTimeSteps = nofLegs * nofTauSteps;
  
%   % object four vector in lab frame (frame [1])
%   traceObj_4v_fr1 = zeros( 4, nofPtrGridObj, nofTimeSteps) * NaN;
  
%   % proper time at each object point
%   timProper = zeros( nofPtrGridObj, nofTimeSteps) * NaN;
    
  % lorentz transform from frame (k) to lab frame (frame (1)),
  % initialize for k=1
  LTToFr1 = eye(4);
  
  
  % number of points used to characterize the event horizon
  % nofHrzPts = 20;
  
  % for each time step (nofTimeSteps) we have  
  %   nofHrzPts (number of points in one horizon) 
  % times
  %   nofLegs (number of boost legs = 5) 
  % 4-vector coordinates
  
  traceHrz_4v_fr1 = zeros( 4, nofHrzPts, nofLegs, nofTimeSteps) * NaN;
    
  % loop over all five boost segments
  % (we do 'nofTauSteps' time steps at once)
  for k = (1:nofLegs),

    disp([mfilename ': leg = ' num2str(k)])
           
%     % Lorentz trafo from current frame to lab frame (frame [1])
%     traceObj_4v_fr1(:,p,tSel) = LTToFr1 * [tmpTim; tmpPos];
    
    tSel = (1:nofTauSteps) + (k-1)*nofTauSteps;
    
%     % four-position of reference point at first time epoch
%     % in frame [k]
%     refPoint = LTToFr1 \ squeeze( traceObj_4v_fr1(:,1,tSel(1)));
    
    % four-position of reference point in frame [k]
    refPoint = LTToFr1 \ squeeze( traceObj_4v_fr1(:,1,tSel));
        
    boostDir = eNrm(:,k);
      
    evHrz4Pos      = zeros( 4, nofTauSteps);

    evHrz4Pos(1,:) = refPoint(1,:);

    % direction perpendicular to boostDir in frame [k]
    boostXDir = [boostDir(2); -boostDir(1); 0];

    % event horizon of current frame
    for p = (1:nofHrzPts),
      
      ofs = boostXDir ...
        * ((p-1) / (nofHrzPts-1) * hrzLength - hrzLength/2);
      
      % use first time epoch
      tmp = refPoint(2:4,1) - boostDir / accRef + ofs;
      
      evHrz4Pos(2:4,:) = repmat( tmp, 1, nofTauSteps);
      
      % transform from frame [k] to frame [1]
      traceHrz_4v_fr1(:,p,k,tSel) = LTToFr1 * evHrz4Pos;
      
    end
        
    % in boost leg k :
    %   copy event horizon position from previous legs (1,...,k-1)
    %   take time from reference proper time
    % *** checkme: is this correct? ***
    % i.e. horizon is stationary in lab frame once next boost is active.
    
    % loop over past ('inactive') event horizons
    tPosFr = (k-2)*nofTauSteps + (1:nofTauSteps);
    tPosTo = (k-1)*nofTauSteps + (1:nofTauSteps);      

    % temporal component is adjusted accordingly
    tTimFr = (k-1)*nofTauSteps + (1:nofTauSteps);
    tTimTo = (k-1)*nofTauSteps + (1:nofTauSteps);      
    
    for kk = (k-1:-1:1)
      
      assert( all( all( all( isnan( traceHrz_4v_fr1(:,:,kk,tPosTo))))) )
      assert( ~any( any( any( isnan( traceHrz_4v_fr1(:,:,kk,tPosFr))))) )
                  
      % spatial part is just copied
      traceHrz_4v_fr1(2:4,:,kk,tPosTo) = traceHrz_4v_fr1(2:4,:,kk,tPosFr);
            
      assert( ~any( any( any( isnan( traceHrz_4v_fr1(:,:,kk+1,tTimFr))))) )
      
      traceHrz_4v_fr1(1,:,kk,tTimTo) = traceHrz_4v_fr1(1,:,kk+1,tTimFr);
      
    end
       
%     % newPosGridObjPts_4v : 4-positions of grid points in next frame
%     posGridObjPts_4v = newPosGridObjPts_4v;
    
%    posEvHrzPts_4v   = newPosEvHrzPts_4v;

    % Lorentz transformation from frame [k] to frame [1] (lab frame)
    % what we want is:
    %   LTToFr1 = LTToFr1 * inv( LT(:,:,k));
    % the following numerically more stable than  
    %   LTToFr1 := LTToFr1 * inv( LT(:,:,k))
    LTToFr1 = LTToFr1 / LT(:,:,k);
    
  end % --- for k = 1:nofLegs, ---

  if 0,
    squeeze( traceHrz_4v_fr1(:,1,:,:))    
    keyboard    
  end
  
  return
        
end

% -------------------------------------------------------------------------
%
%  relrotcosangles5seg  cosines of boost angles for closed loop motion
%  
%  [C12, C23, C13] = relrotcosangles5seg( be)
%
%  closed loop motion : series of boost such that final velocity is zero
%                       and final position is the same as the start point. 
%
%  C12 = cos( alpha12)
%  alpha12 : angle between boost vectors e1 and e2
%
%  Note:
%    alpha12 + alpha23 + alpha13 = 2*pi
%
function [C12, C23, C13, T12, T23] = local_relrotcosangles5seg( be)

  % result for C12 is numerically unstable for small beta
%  assert( all( be > 0.05 & be < 1));

  ga = 1 ./ sqrt( 1 - be.^2);

  [C12, C23, C13] = local_cosinesolution( ga, be);
         
  a3 = ((10*ga + 6) ./ (2*ga + 1));
  a2 = ((2*ga.^3 - 3*ga.^2 + 10*ga + 9) ./ (2*ga + 1));
  a1 = (2*ga .* (ga - 1) .* (4*ga + 3) ./ (2*ga + 1));
  a0 = (- ga.^2 .* (4*ga.^2 + 2*ga - 1) ./ (2*ga + 1));
  
%  z = memo44193( a3, a2, a1, a0);
  T1 = -a3 / 4;
  T2 = a2.^2 - 3*a3.*a1 + 12*a0;
  T3 = (2*a2.^3 - 9*a3.*a2.*a1 + 27*a1.^2 + 27*a3.^2.*a0 - 72*a2.*a0) / 2;
  T4 = (-a3.^3 + 4*a3.*a2 - 8*a1) / 32;
  T5 = (3*a3.^2 - 8*a2) / 48;

  R1 = sqrt( T3.^2 - T2.^3); % power -> 1st primitive root
%  R2 = power( T3 + R1, 1/3);
  R2 = (T3 + R1).^(1/3);
  R3 = (1/12) * (T2 ./ R2 + R2);
%  R4 = power(T5 + R3, 1/2);
  R4 = sqrt(T5 + R3);
  R5 = 2*T5 - R3;
  
  R6 = ones( size( R4));
  idx = (R4 ~= 0);
  R6(idx) = T4(idx) ./ R4(idx);
%   if abs(T4) < 1e-10 && abs(R4) < 1e-10 % ideally: if T4 = R4 = 0
%     R6 = 1;
%     R4 = 0;
%   else
%   end
  
%  z = zeros(1,4);
  
%   z(1) = T1 + R4 + power(R5 + R6, 1/2);
%   z(2) = T1 + R4 - power(R5 + R6, 1/2);
%   z(3) = T1 - R4 + power(R5 - R6, 1/2);
%   z(4) = T1 - R4 - power(R5 - R6, 1/2);
%   T23Sqa = T1 + R4 + sqrt( R5 + R6);
%   T23Sqb = T1 + R4 - sqrt( R5 + R6);
%   T23Sqc = T1 - R4 + sqrt( R5 - R6);
%   T23Sqd = T1 - R4 - sqrt( R5 - R6);
  
  T23Sqa = T1 + sqrt(T5 + R3) + sqrt( R5 + R6);
%  T23Sqb = T1 + sqrt(T5 + R3) - sqrt( R5 + R6);
  T23Sqc = T1 - sqrt(T5 + R3) + sqrt( R5 - R6);
%  T23Sqd = T1 - sqrt(T5 + R3) - sqrt( R5 - R6);
  
  sel = ((T3 + R1) >= 0 );

  T23 = zeros( size( ga)) * NaN;
  
  T23(~sel) = sqrt( T23Sqa(~sel));
  T23(sel)  = sqrt( T23Sqc(sel));

  % max( abs( imag( T23)))
  
  assert( max( abs( imag( T23))) < 1e-4)
    
  T23 = real( T23);
  
  T12 = (-T23.^2 + ga + sqrt(T23.^4 + T23.^2 + ga.^2))./T23;
    
  assert( max( abs( imag( T12))) < 1e-6)

  tmpT12 = T12;
  tmpT23 = T23;
  
  
  [T12, T23] = local_tangenthalfanglesolution( ga);
  
  C12 = ((1-T12.^2) ./ (1+T12.^2));
  C23 = ((1-T23.^2) ./ (1+T23.^2));

%   C12
%   tmpC12
%   C12 - tmpC12
%   C23
%   tmpC23
%   C23 - tmpC23
%   keyboard
  
%   tmpT12
%   T12
%   tmpT12 - T12
%   tmpT23
%   T23
%   tmpT23 - T23
%   
%   keyboard
  
  return
  
  %   clfig(99)
%   plot( be, C12)
%   hold on
%   plot( be, tmpC12, 'r--')
% 
%   clfig(98)
%   plot( be, tmpC12 - C12)
%   keyboard
%   
%   max( abs( tmpC12 - C12))
%   assert( max( abs( tmpC12 - C12)) < 1e-9)

  tmpAA = (ga - 1).*(ga + 1).^3;
  tmpBB = 2*(ga - 1).*(ga + 1).^2 .* (2*ga + 1);
  tmpCC = (ga + 1).*(6*ga.^3 - 2*ga.^2 - ga + 3);
  tmpDD = 2*(ga + 1).*(2*ga.^3 - 3*ga.^2 + 3*ga + 2);
  tmpEE = (  ga.^4 - 2*ga.^3 +   ga.^2 -  4*ga - 4);

  delta = C23.^4 .* tmpAA + C23.^3 .* tmpBB + C23.^2 .* tmpCC + C23 .* tmpDD + tmpEE;
          
  figure(99)
  clf
  plot( be, delta)
  keyboard

  return
  
  
%  C12 = ( C23 .* ga - (1-ga) ./ ((ga.^2-1)./ga.^2) + C13 ) ./ ( (1-ga) .* ( C23 + 1) );
%  C12 = -( C23 .* ga .* (ga+1) + ga.^2 + C13 .* (ga+1) ) ./ ( (ga-1) .* ( C23 + 1) .* (ga+1) );
  tmpC13 = C13;
  tmpC23 = C23;
  
  tmpC13(tmpC13>1)  = 1;
  tmpC13(tmpC13<-1) = -1;
  tmpC23(tmpC23>1)  = 1;
  tmpC23(tmpC23<-1) = -1;
  tmpC12 = tmpC13 .* tmpC23 - sin( acos( tmpC13)) .* sin( acos( tmpC23));
  
%   tmpC12 = ( C23 .* ga - (1-ga) ./ ((ga.^2-1)/ga.^2) + C13 ) ./ ( (1-ga) .* ( C23 + 1) );
% 
%   assert(C12 == tmpC12)
  
%  tmp = ( C23 .* ga .* (ga+1) + ga.^2 + C13 .* (ga+1) );
  figure(99)
  clf
  plot( be, C12)
  hold on
  plot( be, tmpC12, 'r-')

  figure(98)
  clf
  plot( be, tmpC13)
  hold on
  plot( be, tmpC23, 'r-')

  figure(97)
  clf
  plot( be, sqrtArg1)
  hold on
  plot( be, sqrtArg2, 'r-')
  plot( be, sqrtArg3, 'g-')

  figure(96)
  clf
%   plot( be, sqrtArg1)
%   hold on
%   plot( be, sqrtArg2, 'r-')
  plot( be, sqrtArg3, 'g-')
  hold on
  plot( be, tmpSqrtArg3, 'r--')
  
  ylabel('sqrtArg3')
  
  keyboard

end

% ----------------------------------------------------------
%
%
%
%
function C23 = local_check( ga, AA, BB, CC, DD, EE, FF)

  S1 = (ga .* CC) ./ ((ga-1) .* (ga+1));
  
  S4 = (8*sqrt(2) * ga.^2 .* sqrt(S1)) ./ ...
       (3*sqrt(3) * (ga-1).^2 .* (ga+1).^5) - AA ./ BB;
     
  S2 = DD .* S4.^(1/3) + EE + FF ./ S4.^(1/3);
  
  S3 = (12*ga.^2) ./ ((ga+1).*sqrt(S2)) - S4.^(1/3) - ...
    FF ./ (DD .* S4.^(1/3)) + ...
    (2*(4*ga.^2-7*ga-9)) ./ (3*(ga-1).*(ga+1).^2);

  C23 = sqrt( S3) / 2 + sqrt( S2) ./ (6*(ga-1).*(ga+1).^2) - ...
        (2*ga+1) ./ (2*ga+2);
  
end

% -------------------------------------------------------------------------
%
%
%
function [C12, C23, C13] = local_cosinesolution( ga, be)

  a  = ga;

% okay, this is ugly, but it's the correct solution!

%  tmpBB = (27*a.^9 +  81*a.^8 - 216*a.^6 -  162*a.^5 + 162*a.^4 + 216*a.^3 - 81*a - 27);
  BB = 27 * (a - 1).^3 .* (a + 1).^6;
  
%  assert( max( abs( BB - tmpBB)) < 1e-5)
  
  AA = ( 8 * a.^6 + 516 * a.^5 - 678 * a.^4 - 2357 * a.^3 ...
         - 135 * a.^2 + 1701 * a + 729); 
  
%  tmpDD = ( 9*a.^6 +  18*a.^5 -   9*a.^4 -   36*a.^3 -  9*a.^2 + 18*a + 9);
  DD = 9 * (a - 1).^2 .* (a + 1).^4;

%  assert( max( abs( DD - tmpDD)) < 1e-9)
  
%  tmpEE = (12*a.^5 -   9*a.^4 -  60*a.^3 -   18*a.^2 + 48*a + 27);
  EE = 3 * (a - 1) .* (a + 1).^2 .* (4 * a.^2 - 7 * a - 9);

%  assert( max( abs( EE - tmpEE)) < 1e-10)
  
  FF = ( 4 * a.^4 - 116 * a.^3 -  59 * a.^2 +  126 * a + 81);
  
  CC = ( 4 * a.^4 +  28 * a.^3 + 193 * a.^2 +  234 * a + 81);
  
%   C23 = sqrt( (12 .* (a-1) .* a.^2 .* (a+1).^2) ./ ((a.^4+2 .* a.^3-2 .* a-1) ...
%     .* sqrt( ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( (a .* CC) ./ ((a-1) .* (a+1)))) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3) ...
%     .* (DD .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt((a .* CC) ./ ((a-1) .* (a+1)))) ./ ((a-1).^2 .* (a+1).^5) ...
%     - AA ./ BB).^(2 ./ 3) + EE .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( (a .* CC) ./ ((a-1) .* (a+1)))) ./ ((a-1).^2 ...
%     .* (a+1).^5) - AA ./ BB).^(1 ./ 3) + FF))) ...
%     - ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt((a .* CC) ./ ((a-1) .* (a+1)))) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1 ./ 3) ...
%     - (FF .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( (a .* CC) ./ ((a-1) .* (a+1)))) ...
%     ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3)) ./ DD ...
%     - ((6 .* a.^3-2 .* a.^2-a+3) ./ (a.^3+a.^2-a-1)) ./ 3+(2 .* a.^2-5) ./ (a.^2-1)) ./ 2 ...
%     + (sqrt(((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt((a .* CC) ./ ((a-1) .* (a+1)))) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3) ...
%     .* (DD .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) ...
%     .* a.^2 .* sqrt((a .* CC) ./ ((a-1) .* (a+1)))) ...
%     ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(2 ./ 3) ...
%     + EE .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( (a .* CC) ./ ((a-1) .* (a+1)) ) ) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1 ./ 3) + FF)) ...
%     ./ ((a-1) .* (a+1).^2)) ./ 6-(2 .* a+1) ./ (2 .* a+2);

  % note: sqrtArg1 diverges for a -> 1
  sqrtArg1 = (a .* CC) ./ ((a-1) .* (a+1));
  sqrtArg1(sqrtArg1<0) = NaN;

%   C23 = sqrt( (12 .* (a-1) .* a.^2 .* (a+1).^2) ./ ((a.^4+2 .* a.^3-2 .* a-1) ...
%     .* sqrt( ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3) ...
%     .* (DD .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) ...
%     - AA ./ BB).^(2 ./ 3) + EE .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1 ./ 3) + FF))) ...
%     - ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1 ./ 3) ...
%     - (FF .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3)) ./ DD ...
%     - ((6 .* a.^3-2 .* a.^2-a+3) ./ (a.^3+a.^2-a-1)) ./ 3+(2 .* a.^2-5) ./ (a.^2-1)) ./ 2 ...
%     + (sqrt(((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3) ...
%     .* (DD .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ... 
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(2 ./ 3) ...
%     + EE .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1 ./ 3) + FF)) ...
%     ./ ((a-1) .* (a+1).^2)) ./ 6-(2 .* a+1) ./ (2 .* a+2);

%   sqrtArg2 = ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3) ...
%     .* ( DD .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ... 
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(2 ./ 3) ...
%     + EE .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1 ./ 3) + FF);
  
%   tmpSqrtArg2 = ( DD .* ((2^(7/2)*3^(-3/2) .* a.^2 ... 
%          .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1/3) ...
%          + EE ...
%          + FF / (((2^(7/2)*3^(-3/2) .* a.^2 ...
%          .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1/3)));
%   
%   delta = abs( sqrtArg2 - tmpSqrtArg2);
%   clfig(99)
%   plot( delta)
%   keyboard
    
  fct = 2^(7/2) * 3^(-3/2);
  
%   tmpsqrtArg2 = ( (fct .* a.^2 .* sqrt( sqrtArg1)) ./ ...
%                   ((a-1).^2 .* (a+1).^5) ...
%                   - AA ./ BB ).^(-1/3) ...
%       .* (DD .* ( (fct .* a.^2 .* sqrt( sqrtArg1)) ./ ...
%                   ((a-1).^2 .* (a+1).^5) ...
%                   - AA ./ BB).^(2/3) ...
%       + EE .* ( (fct .* a.^2 .* sqrt( sqrtArg1)) ./ ...
%                 ((a-1).^2 .* (a+1).^5) ...
%                   - AA ./ BB).^(1/3) + FF);

  tmp = ( (fct .* a.^2 .* sqrt( sqrtArg1)) ./ ...
                  ((a-1).^2 .* (a+1).^5) ...
                  - AA ./ BB );                
                
%   sqrtArg2 = tmp.^(-1/3) .* ( DD .* tmp.^(2/3) ...
%                               + EE .* tmp.^(1/3) + FF);
  sqrtArg2 = ( DD .* tmp.^(1/3) + EE + FF ./ tmp.^(1/3));
                
                
%   sqrtArg4 = ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^((-1) ./ 3) ...
%     .* (DD .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(2 ./ 3) ...
%     + EE .* ((2.^(7 ./ 2) .* 3.^((-3) ./ 2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1 ./ 3) + FF);
%   
%   assert( all( sqrtArg2 == sqrtArg4))
%   
%   sqrtArg4(sqrtArg4<0) = NaN;

  sqrtArg2(sqrtArg2<0) = NaN;

  sqrtArg4 = sqrtArg2;
  
%   tmpSqrtArg3 = (12 .* (a-1) .* a.^2 .* (a+1).^2) ...
%     ./ ((a.^4 + 2 .* a.^3 - 2 .* a - 1) .* sqrt( sqrtArg4) ) ...
%     - ((2.^(7/2) .* 3.^(-3/2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1/3) ...
%     - (FF .* ((2.^(7/2) .* 3.^(-3/2) .* a.^2 ...
%     .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(-1/3)) ./ DD ...
%     - ((6 .* a.^3 - 2 .* a.^2 - a + 3) ./ (a.^3 + a.^2 - a - 1)) ./ 3 + (2 .* a.^2-5) ./ (a.^2-1);

  % 2^(7/2)*3^(-3/2) = sqrt(128/27)
%  fct = 2^(7/2)*3^(-3/2);

%   sqrtArg3 = (12 .* a.^2 ) ./ ( (a+1) .* sqrt( sqrtArg4) ) ...
%       - ( (fct .* a.^2 .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(1/3) ...
%     - ( FF .* ...
%         ( (fct .* a.^2 .* sqrt( sqrtArg1)) ./ ((a-1).^2 .* (a+1).^5) - AA ./ BB).^(-1/3)) ./ DD ...
%     - ( (6 .* a.^3 - 2 .* a.^2 - a + 3) ./ ((a-1) .* (a+1).^2)) ./ 3 + (2 .* a.^2 - 5) ./ (a.^2 - 1);
  
  sqrtArg3 = (12 .* a.^2 ) ./ ( (a+1) .* sqrt( sqrtArg4) ) ...
      - (sqrtArg2 - EE) ./ DD ...  
      - ( (6 .* a.^3 - 2 .* a.^2 - a + 3) ...
      ./ ((a-1) .* (a+1).^2)) ./ 3 + (2 .* a.^2 - 5) ./ (a.^2 - 1);
  
  sqrtArg3(sqrtArg3<0) = NaN;
  
  % note: the following is undefined for a -> 1
  C23 = sqrt( sqrtArg3 ) ./ 2 ...
        + ( sqrt( sqrtArg2) ./ ((a-1) .* (a+1).^2)) ./ 6 ...
        - (2*a + 1) ./ (2*a + 2);

% check formulas of paper       
%   tmpC23 = local_check( ga, AA, BB, CC, DD, EE, FF);  
%   delta = tmpC23 - C23;
%   
%   clfig(99)
%   plot( be, delta)
%   ylabel('delta')
%   xlabel('\beta')
%   keyboard
            
  C13 = ( -ga ./ (ga+1) ...
          .* (C23.^3 .* ga.^2 - C23.^3 + 2*C23.^2 .* ga.^2 - 2*C23.^2 ...
              + C23 .* ga.^2 + C23 .* ga + 2*C23 + ga + 1) ...
          ./ (C23.^2 .* ga.^2 - C23.^2 + 2*C23 .* ga.^2 - C23 .* ga ...
            - C23 + ga.^2 - ga + 2));

  
%  tmpC12 = ( C23 .* ga - (1-ga) ./ be.^2 + C13 ) ./ ( (1-ga) .* ( C23 + 1) );
  C12 = -( C23 .* ga .* (ga+1) + ga.^2 + C13 .* (ga+1) ) ...
          ./ ( (ga-1) .* ( C23 + 1) .* (ga+1) );

        
  sel = find( be == 0);
  
  C12(sel) = -17/18 + sqrt(10)/18;
  
  C23(sel) = (-11/4*(-14*sqrt(10)+247) ...
    - 1/2*sqrt(-140*sqrt(10)+2470) ...
    + 107*sqrt(5)/12 + 1165*sqrt(2)/12) ...
    / (sqrt(-14*sqrt(10)+247)*(2*sqrt(10)+11));
  
  C13(sel) = NaN;
  
end

% -------------------------------------------------------------------------
%
%  calculate angles between boost directions e1<->e2 and e2<->e3
%
function [T12, T23] = local_tangenthalfanglesolution( ga)

  a = 1;
  b = (4*ga + 8);
  
%   c = (-8*ga.^3 - 36*ga.^2 - 52*ga - 18);
%   d = (-32*ga.^4 - 112*ga.^3 - 104*ga.^2 - 20*ga + 8);
%   e = - 32*ga.^5 - 64*ga.^4 - 40*ga.^3 - 4*ga.^2 + 4*ga + 1;
  
%   Delta = -524288 * ga ...
%                  .* (ga - 1).^3 ...
%                  .* (ga + 1).^7 ...
%                  .* (4*ga.^4 + 28*ga.^3 + 193*ga.^2 + 234*ga + 81);
               
  Delta0 = 16   * (ga + 1).^2 ...
                .* (4*ga.^4 + 28*ga.^3 + 157*ga.^2 + 126*ga + 9);
              
%   Delta1 = -128 * (ga + 1).^3 ...
%                 .* (8*ga.^6 + 84*ga.^5 + 618*ga.^4 + 2071*ga.^3 ...
%                     + 2349*ga.^2 + 729*ga - 27);
                   
  p =  -2 * (ga + 1)    * (4*ga.^2 + 17*ga + 21);
  q = -16 * (ga + 1).^2 * (ga.^2 - 2*ga - 9);
  
  Q = 4 * (-8*ga.^9 - 108*ga.^8 - 894*ga.^7 - 4185*ga.^6 - 10500*ga.^5 ... 
           - 14607*ga.^4 - 11278*ga.^3 - 4455*ga.^2 - 648*ga ...
           + 12 * sqrt(6) ...
                *sqrt(4*ga.^15 + 44*ga.^14 + 317*ga.^13 + 1058*ga.^12 ...
                      + 1316*ga.^11 - 910*ga.^10 - 4275*ga.^9 - 3500*ga.^8 ... 
                      + 1780*ga.^7 + 4720*ga.^6 + 2311*ga.^5 - 854*ga.^4 ...
                      - 1372*ga.^3 - 558*ga.^2 - 81*ga) ...
           + 27).^(1/3);

  S = (1/2) * sqrt( - (2/3) .* p + 1 ./ (3*a) .* (Q + Delta0 ./ Q));
  
%   T12Sq = (- b ./ (4*a) - S + (1/2) * sqrt(-4*S .^ 2 - 2 * p + q ./ S));
%   T12Sq = (- b ./ (4*a) - S - (1/2) * sqrt(-4*S .^ 2 - 2 * p + q ./ S));
  % >>> this is the correct solution! <<<
  T12Sq = (- b ./ (4*a) + S + (1/2) * sqrt(-4*S .^ 2 - 2 * p - q ./ S));
%  T12Sq = (- b ./ (4*a) + S - (1/2) * sqrt(-4*S .^ 2 - 2 * p - q ./ S));

%   T12Sq
%   tmp1 = T12Sq.^4 .* a + T12Sq.^3 .* b + T12Sq.^2 .* c + T12Sq .* d + e 

%   % alternative real solution
%   T12Sq = (- b ./ (4*a) + S - (1/2) * sqrt(-4*S .^ 2 - 2 * p - q ./ S));

%   T12Sq
%   tmp2 = T12Sq.^4 .* a + T12Sq.^3 .* b + T12Sq.^2 .* c + T12Sq .* d + e
%   
%   keyboard
  
%   x1 = (- b./(4*a) - S + (1/2) * sqrt(-4*S.^2 - 2*p + q./S));
%   x2 = (- b./(4*a) - S - (1/2) * sqrt(-4*S.^2 - 2*p + q./S));
%   x3 = (- b./(4*a) + S + (1/2) * sqrt(-4*S.^2 - 2*p - q./S));
%   x4 = (- b./(4*a) + S - (1/2) * sqrt(-4*S.^2 - 2*p - q./S));

  T12 = sqrt( T12Sq);

  assert( max( abs( imag( T12))) < 1e-8);
  T12 = real( T12);
  
  T23 = ((-2*T12.*(ga + 1) ... 
       + sqrt(-T12.^4 + 8*T12.^2*ga + 6*T12.^2 + 8*ga.^3 + 8*ga.^2 - 1)) ...
       ./ (-T12.^2 + 2*ga + 1));

%   % second solution   
%  T12 = +sqrt( 4*ga + 2*sqrt(2).*sqrt((ga + 1).^3) + 3);
%  T12 = +sqrt( 4*ga - 2*sqrt(2).*sqrt((ga + 1).^3) + 3);
%  T12 =+sqrt( 4*ga + 2*sqrt(2).*sqrt((ga + 1).^3) + 3);
%  T12 = -sqrt( 4*ga + 2*sqrt(2).*sqrt((ga + 1).^3) + 3);
%   T12a = +sqrt( 4*ga - 2*sqrt(2)*sqrt((ga + 1)**3) + 3)
%   T12a = -sqrt( 4*ga + 2*sqrt(2)*sqrt((ga + 1)**3) + 3)
%   T12a = -sqrt( 4*ga - 2*sqrt(2)*sqrt((ga + 1)**3) + 3)
%   % second solution   
%  assert( max( abs( imag( T12))) < 1e-8);
%  T12 = real( T12);
%  
%  T23 = ((-2*T12.*(ga + 1) ... 
%       - sqrt(-T12.^4 + 8*T12.^2*ga + 6*T12.^2 + 8*ga.^3 + 8*ga.^2 - 1)) ...
%       ./ (-T12.^2 + 2*ga + 1));

  assert( max( abs( imag( T23))) < 1e-8);
  T23 = real( T23);
  
end

% -------------------------------------------------------------------------
%
% ISODD True for odd matrix elements.
%
%   y = isodd( x)
%
%	returns 1 for odd elements of X  and 0 otherwise.
%
function y = local_isodd( x)

  y = ((round((x+1)/2) == (x+1)/2) & isfinite(x));

end
