%
% This computes the effective mass and position of the 2-soliton
% with external potential b(x,t) 
% The input is a function fun=fun(x,t), the time vector T, and initial 
% value ZZ=[a1, a2, c1,c2];
% It returns T, the time vector (which is the same as the time
% input vector t), and Y = [Y(:,1),Y(:,2),Y(:,3),Y(:,4)] which 
% corresponds to [ a1 , a2, c1, c2] at times T(j).
% 
% If the last argument is set to 1 (the default value) then a waterfall plot
% is produced.
%
% Here is an example of an input
%
%   [T,Y]=effdyn(@(x,t) 50*cos(x), linspace(0,0.025,20), [-1,-2,6,9]) 
%

function [T, Y] = effdyn(fun,TT,ZZ,del)

  options = odeset('RelTol',1e-5,'AbsTol',[1e-5 1e-5 1e-5 1e-5]);
  [T,Z] = ode45(@ff,TT,ZZ,options,fun);
  Y = real(Z);

  if nargin < 4, del = 1; end
  if del ~= 1,   return;  end

  x = linspace(-pi,pi,500);
  hold on 
  y1=0;
  y2=0;
  for j=1:length(T)
    q = qu(x,Y(j,:));
    plot3(x,T(j)+0*x,q,'r','LineWidth',1);
    y1=min(min(q),y1);
    y2=max(max(q),y2);
  end
  xlabel x, ylabel t, axis([-pi pi min(T) max(T) y1 y2 ]), grid off
  view(0,60)
  hold off

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

function dy = ff(t,y,fun)

  a1 = sign(y(3));
  a2 = sign(y(4));
  dy = [y(3).^2 - a1*diffB2(t,y,3,fun);
        y(4).^2 - a2*diffB2(t,y,4,fun);
        a1*diffB2(t,y,1,fun);
        a2*diffB2(t,y,2,fun)];
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function qu = qu(x,ZZ)

  k1 = ZZ(3)*(x - ZZ(1));
  k2 = ZZ(4)*(x - ZZ(2));

  ga1 =  exp(-k1);
  ga2 = -exp(-k2);

  m11 = (1+ga1.^2)./(2*ZZ(3));
  m22 = (1+ga2.^2)./(2*ZZ(4));
  m12 = (1+ga1.*ga2)./(ZZ(3)+ZZ(4));

  M  = m11.*m22 - m12.^2;
  M1 = ga1.*(m12 - m22) + ga2.*(m12-m11);
  qu = M1./M;

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

function dB2 = diffB2(t,W,j,fun)
  
  h = 1e-5;
  Wp = W;  Wp(j) = Wp(j) + h;
  Wm = W;  Wm(j) = Wm(j) - h;
  
  x   = linspace(-6,6,4000);
  up  = qu(x,Wp);
  um  = qu(x,Wm);
  y   = fun(x,t).*(up-um).*(up+um)/2;
  dB2 = trapz(x,y)/(2*h);

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