%   Movie Project for EPS 109;  Version 2.  (Fixed energy calcs.)
%   Nathan Haese, November 2009
%
%   The project is to follow the motion of two positive ions in the 
%   field of a dust grain with high negative charge.
%   First ion is carbon ion.  Second ion is protonated hydrogen, H3+.

    clear; clf;
    global e1 qd q1 q2 md m1 m2
    
%   Physical constants
    Navo    =   6.0221415e23;           % Avogadro's number
    kboltz  =   1.380654e-23;           % Boltzmann constant J/deg
    e0      =   8.854187e-12;           % Permittivity in F/m
    e1      =   1/4/pi/e0;              % convenient constant in m/F
    echg    =   1.602176487e-19;        % mag of electron charge in coulombs
%   Gas conditions                      
    T       =   5;                     % temperature in K
%   Dust particle specifications
    nqd    =   1000;                   % no. of electron chg on dust
    qd     =  -nqd*echg;               % charge in coulombs
    rd     =   1e-07;                  % radius in m
    rhod   =   5;                      % density of dust gram/cc
    rhod   =   rhod*1000;              % density in Kg/m^3
    md     =   rhod*(4*pi/3)*rd^3;     % mass of dust particle in Kg
    xd     =   0;                      % x-coord of dust particle
    yd     =   0;                      % y-coord of dust particle
    vxd    =   0;                      % vx of dust particle
    vyd    =   0;                      % vy of dust particle
%   Ion_1 specifications                
    mw1     =   12;                      % molecular weight grams/mole
    nq1     =   1;                      % no. of electron chg on ion
    m1      =   mw1/Navo/1e03;          % mass in Kg
    q1      =   nq1*echg;               % charge in coulombs   
    r1      =   1776e-06;                      % polar dust-ion dist. in meters
    phi1    =   0;                      % polar angle to ion in deg
    x1      =   r1*cos(phi1*pi/180);    % x-coord of ion
    y1      =   r1*sin(phi1*pi/180);    % x-coord of ion
    theta1  =   90;                     % pol. angle of ion veloc vec (deg)
    v1      =   sqrt(3*kboltz*T/m1);    % ion velocity in m/s
    vx1     =   v1*cos(theta1*pi/180);  % vx of ion-1
    vy1     =   v1*sin(theta1*pi/180);  % vy of ion-1
%   Ion-2 specifications                
    mw2     =   3;                      % molecular weight grams/mole
    nq2     =   1;                      % no. of electron chg on ion
    m2      =   mw2/Navo/1e03;          % mass in Kg
    q2      =   nq2*echg;               % charge in coulombs   
    r2      =   1*r1;                      % polar dust-ion dist. in meters
    phi2    =   180;                    % polar angle to ion in deg
    x2      =   r2*cos(phi2*pi/180);    % x-coord of ion
    y2      =   r2*sin(phi2*pi/180);    % x-coord of ion
    theta2  =   270;                          % pol. angle of ion veloc vec (deg)
    v2      =   sqrt(3*kboltz*T/m2);     % ion velocity in m/s
    vx2     =   v2*cos(theta2*pi/180);  % vx of ion
    vy2     =   v2*sin(theta2*pi/180);  % vx of ion 
%   Analytical Orbit Parameters of single ion   % From Newton's equations 
    vc1      =   sqrt(abs(e1*qd*q1)/m1/r1);   % critical velocity
    vc2      =   sqrt(abs(e1*qd*q2)/m2/r2);   % critical velocity
    r0      =   1e06*abs(e1*qd*q1/3/kboltz/T);     % critical rad in mu
    eccen1   =   (v1/vc1)^2 - 1;              % eccentricity
    eccen2   =   (v2/vc2)^2 - 1;              % eccentricity
                %    plot(w(1,1),w(2,1),'or')
                %    phi2 = (pi/180)*(0:1:359);
                %    re = r_ion * (v_ion/vc)^2./(1 + ((v_ion/vc)^2-1).*cos(phi2));
                %    plot(re.*cos(phi2), re.*sin(phi2),'-r')
                %    axis square
%   Calculate Center-of-Mass    
    xcom = (xd*md + x1*m1 + x2*m2)/(md+m1+m2);
    ycom = (yd*md + y1*m1 + y2*m2)/(md+m1+m2);
    
% initialize array of variables for RK calculation
    w(1,1)  =   xd - xcom;      
    w(2,1)  =   yd - ycom;      
    w(3,1)  =   vxd;
    w(4,1)  =   vyd;
    w(5,1)  =   x1 - xcom;
    w(6,1)  =   y1 - ycom; 
    w(7,1)  =   vx1;
    w(8,1)  =   vy1;
    w(9,1)  =   x2 - xcom;
    w(10,1)  =  y2 - ycom; 
    w(11,1)  =  vx2;
    w(12,1)  =  vy2;
    
% integration
    [t,w] = ode45('wdotHaese2', [0,0.005], w); % fast, variable time step ode function

% creat arrays for energies and angular momentum 
    N = length(w);     % number of steps used by ode23
    K = zeros(N,1);    % kinetic energy
    U = zeros(N,1);    % potential energy
    E = zeros(N,1);    % total energy
    L = zeros(N,1);    % angular momentum
 
% movie commands 
    dum = datestr(now); % dummy variable for naming files
    movie_file = ['Haese_', dum(13:14),dum(16:17),dum(19:20)];
    set(gca, 'nextplot','replacechildren');
    aviobj = avifile(movie_file,'quality',100,'compression','Cinepak');

 % Now, compute the energy quantities.
    for i = 1:N
        K(i,1)  = 0.5*(md*(w(i,3)^2 + w(i,4)^2) + m1*(w(i,7)^2+w(i,8)^2) + m2*(w(i,11)^2+w(i,12)^2));
        r2d1    = sqrt(((w(i,1)-w(i,5))^2 + (w(i,2)-w(i,6))^2));
        r2d2    = sqrt(((w(i,1)-w(i,9))^2 + (w(i,2)-w(i,10))^2));
        r212    = sqrt(((w(i,5)-w(i,9))^2 + (w(i,6)-w(i,10))^2));
        U(i,1)  = e1*qd*q1/r2d1 + e1*qd*q2/r2d2 + e1*q1*q2/r212;
        E(i,1)  = K(i) + U(i);
        L(i,1)  = md*(w(i,1)*w(i,4)-w(i,2)*w(i,3))+...
               m1*(w(i,5)*w(i,8)-w(i,6)*w(i,7))+...
               m2*(w(i,9)*w(i,12)-w(i,10)*w(i,11));
    end
     
   % subplot(2,1,1)
   % plot(t,K,t,U,t,E,t,L), legend('K','U','E','L'), xlabel('t(s)')
   % subplot(2,1,2) 
    
    % coordinate rotation to dust to ion-1 frame (not in movie)
         
    for i=1:length(w)
        phi3    = atan2(w(i,6), w(i,5));  %rotation angle
        c       = cos(phi3);  % rotation factor
        s       = sin(phi3);  % rotation factor
        xdt     = c*w(i,1) + s*w(i,2);       
        ydt     = -s*w(i,1) + c*w(i,2);
        x1t     = c*w(i,5) + s*w(i,6);       
        y1t     = -s*w(i,5) + c*w(i,6);
        x2t     = c*w(i,9) + s*w(i,10);       
        y2t     = -s*w(i,10) + c*w(i,10);
     
        plot(w(i,1), w(i,2),'oR', w(i,5), w(i,6),'o', w(i,9), w(i,10),'o'), axis ([-0.02 0.02 -0.02 0.02])
        xlabel('x(m)'), ylabel('y(m)')
        F = getframe(gca);
        aviobj = addframe(aviobj, F);
        pause(0.01)
                %  plot(xdt, ydt, 'o', x1t, y1t,':', x2t, y2t,'-'), axis ([-0.002 0.002 -0.002 0.002]), hold on
                %  pause(0.03)
    end
    xlabel('x(m)'), ylabel('y(m)')
    aviobj = close(aviobj);
    
 

    
    
    
        
    
    