function simulation_run( particle_startVel, particle_startPos, handles )
% simulation loop



%     Copyright (C) 2007 
%
%       mag. David Erzen
%       Faculty of Mechanical Engineering
%       LECAD Laboratory
%       Askerceva 6
%       1000 Ljubljana
%       SLOVENIA
%       contact email: david.erzen@lecad.uni-lj.si
%       
%       Prof. John P. Verboncoeur
%       Plasma Theory and Simulation Group
%       University of California
%       Berkeley, CA 94720-1730 USA
%       
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
    


    global_variables;
    
    status = STATUS_RUNNING;            % simulation_running is set to 0 by user to stop or pause the simualtion
    particle_vel = particle_startVel;    % temporary velocity and position information for solver 
    particle_pos = particle_startPos;
    trajectory(1, 1:3) = particle_startPos;
    w = vrworld('sim1.wrl');
    
    counterSim = 0;
    absolute_time = 0;
    
    % SIMULATION LOOP
    while status == STATUS_RUNNING
    
       % calcuate to (tspan) for solver 
       % TODO: to should be calculated as 1/100 part of inverse value of
       % Deby frequency
        to = time_scale( particle_pos );        
        absolute_time = absolute_time + to;
        set(w, 'Time', absolute_time*1000);
        
        % calcualte one step with function field_Efunction and
        % field_Bfunction
        % TODO: solver calculates more points then needed afterwards, this
        % should not be
        
        [ time, partial_trajectory ] = ode15s( @acceleration, [0, to], [particle_pos, particle_vel] );
        
        % update particle_vel and particle_vel
        particle_pos = partial_trajectory ( length( partial_trajectory ), 1:3 );
        particle_vel = partial_trajectory ( length( partial_trajectory ), 4:6 );
        
        
        % update trajectory and guiding center        
        trajectory = cat( 1, trajectory, particle_pos ); 
        calculate_guideCenter;
                
        % update virtual world:
        
        % update force vector
        if vectorForce_set
            % calculate force vector if shown
            acc = acceleration( 0, partial_trajectory ( length( partial_trajectory ), : ) );
            vector_force = (acc(4:6))'/factor_force;
            [ vector_forceRotation, vector_forceConePos, vector_forcePos, vector_forceHeight] = vectorParameters( particle_pos, vector_force );
            w.a_cone_Material.transparency = 0;
            w.a_cyl_Material.transparency = 0;
            w.a_cylinder.rotation = vector_forceRotation;
            w.a_cylinder.translation = vector_forcePos;
            w.a_Cyl.height = vector_forceHeight;
            w.a_cone.rotation = vector_forceRotation;
            w.a_cone.translation = vector_forceConePos;
        else
            % otherwise hide it by putting tranparency 100
            w.a_cone_Material.transparency = 100;
            w.a_cyl_Material.transparency = 100;
        end;

        % update velocity vector
        if vectorVel_set
            % calculate velocity vector
            vector_velocity = particle_vel/factor_velocity;
            [ vector_velocityRotation, vector_velocityConePos, vector_velocityPos, vector_velocityHeight] = vectorParameters( particle_pos, vector_velocity );
            w.v_cone_Material.transparency = 0;
            w.v_cyl_Material.transparency = 0;
            w.v_cylinder.rotation = vector_velocityRotation;
            w.v_cylinder.translation = vector_velocityPos;
            w.v_Cyl.height = vector_velocityHeight;
            w.v_cone.rotation = vector_velocityRotation;
            w.v_cone.translation = vector_velocityConePos;
        else
            % otherwise hide it by putting tranparency 100
            w.v_cone_Material.transparency = 100;
            w.v_cyl_Material.transparency = 100;
        end;

        %update guiding center vector
        if (vectorGuide_set) & length(guide_center) > 0
            % calculate guiding center vector
            vector_guidingCenter = guide_center(length(guide_center(:,1)), 1:3) - particle_pos;                
            [ vector_guideRotation, vector_guideConePos, vector_guidePos, vector_guideHeight] = vectorParameters( particle_pos, vector_guidingCenter );
            w.g_cone_Material.transparency = 0;
            w.g_cyl_Material.transparency = 0;
            w.g_cylinder.rotation = vector_guideRotation;
            w.g_cylinder.translation = vector_guidePos;
            w.g_Cyl.height = vector_guideHeight - w.g_Cone.height;
            w.g_cone.rotation = vector_guideRotation;
            % adjusting the height of the cylinder and the position of the
            % cone in order to place the tip of the cone in the guiding
            % center
            if vector_guideConePos ~= vector_guidePos
                w.g_cone.translation = vector_guideConePos - w.g_Cone.height*(vector_guideConePos - vector_guidePos)/(2*norm(vector_guideConePos - vector_guidePos));
            else
                w.g_cone.translation = vector_guideConePos;
            end;
        else
            % otherwise hide it by putting tranparency 100
            w.g_cone_Material.transparency = 100;
            w.g_cyl_Material.transparency = 100;
        end;

        % update trajectory
        if traj_set
            w.trajectory_spine.spine = trajectory;
        else
            w.trajectory_spine.spine = trajectory( length( trajectory ), 1:3 );
        end;

        % update guiding center
        if guide_set
            w.guide_spine.spine = guide_center;
        else
            if length( guide_center ) > 0
                w.guide_spine.spine = guide_center( 1, 1:3 );
            end;
        end;

        % update particle
        w.ball.translation = particle_pos;

        % update magnetic field lines for time factor
        switch get( handles.popMagField,'Value' )
            case BFIELD_TIMECONSTANT
            otherwise
                % TODO: get rid of these 18 and 9 asap!!!!!!!!!!!!!!!!!!!!
                for i = 1:18              % updating magnetic field lines    
                    if i <= 9 
                        tmp_num = i;
                    else
                        tmp_num = i+1;
                    end
                    name_line = strcat('B', num2str(tmp_num), '_Material');
                    line = vrnode(w, name_line); 
                    t = absolute_time;
                    line.transparency = 1-abs(eval(field_BTimefunctionX)+eval(field_BTimefunctionY)+eval(field_BTimefunctionZ))/3;                    
                end
        end
         
        % update electric field lines for time factor
        switch get( handles.popElecField,'Value' )
            case EFIELD_TIMECONSTANT
            otherwise
                % TODO: get rid of these 18 and 9 asap!!!!!!!!!!!!!!!!!!!!
                for i = 1:18              % updating electric field lines    
                    if i <= 9
                        tmp_num = i;
                    else
                        tmp_num = i+1;
                    end
                    name_line = strcat('E', num2str(tmp_num), '_Material');
                    line = vrnode(w, name_line); 
                    t = absolute_time;
                    line.transparency = 1-abs(eval(field_ETimefunctionX)+eval(field_ETimefunctionX)+eval(field_ETimefunctionX))/3;    
                end
        end
            
            
            
            vrdrawnow;
        
        
    end
end

% ******************* calculates time scale for one iteration *************
function to = time_scale( pos )

 %   global_variables;
    
%    B = norm( calculate_B( 0, pos ) );
%    if B > 0
%        to = particle_mass/( particle_charge*B );        
%    else
        to = 1e-7;
%    end;    
% TODO: calculate to taking in account the magnetic and electric field
end

% *************** calculates guiding center *******************************
function calculate_guideCenter
% TODO: check if the calculation for the guide center is correct
    global_variables;    
   
    traj_length = length(trajectory);
    if (traj_length > 0 ) traj_length = length(trajectory(:,1)); 
    end;
    if (traj_length > 2)  % trajectory is longer than three points
        t = trajectory(traj_length,1:3);
        t1 = trajectory(traj_length - 1,1:3);
        t2 = trajectory(traj_length - 2,1:3);
        
    if ((t(1) ~= t1(1)) || (t(2) ~= t1(2)) || (t(3) ~= t1(3))) && ((t(1) ~= t2(1)) || (t(2) ~= t2(2)) || (t(3) ~= t2(3))) && ((t2(1) ~= t1(1)) || (t2(2) ~= t1(2)) || (t2(3) ~= t1(3))) % the last three trajectory points are different from one to another            
        s12 = norm(trajectory(traj_length - 1,1:3)-trajectory(traj_length - 2,1:3));
            s23 = norm(trajectory(traj_length ,1:3)-trajectory(traj_length - 1,1:3));
            s13 = s12 + s23;
            
            % calculate vector for geometric guide center vector direction
            temp_guide = 2*(s12.*(trajectory(traj_length ,1:3)-trajectory(traj_length - 1,1:3)) - s23.*(trajectory(traj_length - 1 ,1:3)-trajectory(traj_length - 2,1:3)))/(s13*s23*s12);
            
            % calculate vector for guiding center from origin
            if norm(temp_guide) ~= 0           
                % calculate length of a guide center vector
                mag_field = calculate_B(0, trajectory(traj_length,1:3));
                if norm(mag_field) > 0
                    % calculate perpendicular velocity and length of the
                    % larmour radious
                    perpendicular_vel = particle_vel - mag_field'*(particle_vel*mag_field)/((norm(mag_field))^2);
                    guide_length = abs(norm(perpendicular_vel)/(norm(mag_field)*particle_charge/particle_mass));
                else
                    guide_length = 0;
                end;
                %**********************************************************               
                % larmour radious correction
                temp_guide1 = trajectory(traj_length -1, 1:3) + temp_guide*guide_length/(norm(temp_guide));
                %**********************************************************                
                % geometricly defined center of temporary circular motion
                %temp_guide1 = trajectory(traj_length -1, 1:3) + temp_guide/(norm(temp_guide))^2;
                %**********************************************************                
            else
                temp_guide1 = trajectory(traj_length -1, 1:3);
            end;

            if length(guide_center) == 0
                guide_center(1,1:3) = temp_guide1(1:3);
            else
                guide_center(length(guide_center(:,1))+1, 1:3) = temp_guide1(1:3);   
            end;  
        end;
    end;
end

% ******calculate parameters for force and velocity vector display ********
function [ vector_rotation, vector_conePos, vector_pos, vector_height] = vectorParameters( particle_pos, vector )
    
    if norm( vector ) > 0
        vector_normalised = vector/norm( vector );
    else
        vector_normalised = vector;
    end;
  
    rotation_angle = acos( dot( vector_normalised, [ 0; 1; 0 ] ) );
    rotation_axis = cross( [ 0 1 0 ], vector_normalised );
    vector_rotation = cat( 2, rotation_axis, rotation_angle );
    
    vector_conePos = particle_pos + vector;
    vector_pos = particle_pos + vector/2;
    vector_height = norm( vector );    
end