how to convert matlab file in octave file

[Carletto Rossi]

Hi guys. I would to convert this Matlab script in an Octave file:

N = 4; %number of particles
n = 1000; %size of cell
iteration = 500; %times of measurement
step = 20; %number of steps per-iteration
T = 1.0; %temperature
K = 1.0; %Bolztmann constant
x = zeros(N,1);
y = zeros(N,1);
p = 1;
too_near = 0;
% ###################### initiation ############################
while(p <= N)
   %generate trial position
   x_trial = rand*n;
   y_trial = rand*n;
   for q = 1:p
       d = sqrt( (x_trial-x(q))^2 + (y_trial-y(q))^2 );
       if(d <= 1.0) %the new position is too "near"
           too_near = 1;
           break; %there's no need to continue the loop
       else
           too_near = 0;
       end

    end
    if(too_near == 0) %the new position is "far", accept it
        x(p) = x_trial;
        y(p) = y_trial;
        p = p+1; %step increment
    end
end
% ######################### metropolis #############################
for k=1:iteration
    W_initial = calculate_potential(N,x,y);
    for j=1:step %how many steps per-iteration
        for i=1:N     %for all particles
           %step range between [-0.5,0.5]
           trial_x = rand-0.5;
           trial_y = rand-0.5;
           %save the old values lest the step is not accepted
           x_temp = x(i);
           y_temp = y(i);
           %trial step, bounded by the periodic boundary condition
           x(i) = x(i) + trial_x;
           if (x(i) > n)
               x(i) = x(i) - n;
           elseif (x(i) < 0.0)
               x(i) = n + x(i);
           end
           y(i) = y(i) + trial_y;
           if (y(i) > n)
               y(i) = y(i) - n;
           elseif (y(i) < 0.0)
               y(i) = n + y(i);
           end

            W_final = calculate_potential(N,x,y);
            if( (W_final-W_initial) < 0.0 ) %move is accepted
                W_initial = W_final;
            else %apply metropolis criterion
                if (rand < exp(W_initial-W_final)/(K*T))
                    W_initial = W_final;
                else %move is not accepted, return to initial position
                    x(i) = x_temp;
                    y(i) = y_temp;
                end
            end
        end
    end
    scatter(x,y,3,[0 0 1], 'filled') %plot 2D scatter-graph
    F(k) = getframe; %capture every plot
end
save simulation.mat F;
movie2avi(F, 'simulation.avi');
% ################# calculate potential ###################
function E = calculate_potential(N,x,y)

    E = 0;

     for i=1:N %for all particles
         for j=(i+1):N %compute the potential of j with respect to i
             d = sqrt( (x(j)-x(i))^2 + (y(j)-y(i))^2);
             U = d^(-12) - d^(-6); %Lennard-Jones potential
             E = E + U;
         end
     end

Thanks in advance


_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[James Sherman Jr.-2]

Not to sound obvious, but have you tried running this script in Octave?  If so, what errors do you get?

On Fri, Jul 3, 2009 at 7:14 PM, Carletto Rossi <nuovodna@...> wrote:

Hi guys. I would to convert this Matlab script in an Octave file:

N = 4; %number of particles
n = 1000; %size of cell
iteration = 500; %times of measurement
step = 20; %number of steps per-iteration
T = 1.0; %temperature
K = 1.0; %Bolztmann constant
x = zeros(N,1);
y = zeros(N,1);
p = 1;
too_near = 0;
% ###################### initiation ############################
while(p <= N)
   %generate trial position
   x_trial = rand*n;
   y_trial = rand*n;
   for q = 1:p
       d = sqrt( (x_trial-x(q))^2 + (y_trial-y(q))^2 );
       if(d <= 1.0) %the new position is too "near"
           too_near = 1;
           break; %there's no need to continue the loop
       else
           too_near = 0;
       end

    end
    if(too_near == 0) %the new position is "far", accept it
        x(p) = x_trial;
        y(p) = y_trial;
        p = p+1; %step increment
    end
end
% ######################### metropolis #############################
for k=1:iteration
    W_initial = calculate_potential(N,x,y);
    for j=1:step %how many steps per-iteration
        for i=1:N     %for all particles
           %step range between [-0.5,0.5]
           trial_x = rand-0.5;
           trial_y = rand-0.5;
           %save the old values lest the step is not accepted
           x_temp = x(i);
           y_temp = y(i);
           %trial step, bounded by the periodic boundary condition
           x(i) = x(i) + trial_x;
           if (x(i) > n)
               x(i) = x(i) - n;
           elseif (x(i) < 0.0)
               x(i) = n + x(i);
           end
           y(i) = y(i) + trial_y;
           if (y(i) > n)
               y(i) = y(i) - n;
           elseif (y(i) < 0.0)
               y(i) = n + y(i);
           end

            W_final = calculate_potential(N,x,y);
            if( (W_final-W_initial) < 0.0 ) %move is accepted
                W_initial = W_final;
            else %apply metropolis criterion
                if (rand < exp(W_initial-W_final)/(K*T))
                    W_initial = W_final;
                else %move is not accepted, return to initial position
                    x(i) = x_temp;
                    y(i) = y_temp;
                end
            end
        end
    end
    scatter(x,y,3,[0 0 1], 'filled') %plot 2D scatter-graph
    F(k) = getframe; %capture every plot
end
save simulation.mat F;
movie2avi(F, 'simulation.avi');
% ################# calculate potential ###################
function E = calculate_potential(N,x,y)

    E = 0;

     for i=1:N %for all particles
         for j=(i+1):N %compute the potential of j with respect to i
             d = sqrt( (x(j)-x(i))^2 + (y(j)-y(i))^2);
             U = d^(-12) - d^(-6); %Lennard-Jones potential
             E = E + U;
         end
     end

Thanks in advance


_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave