2 views (last 30 days)
Amna on 5 Dec 2013
Edited: sixwwwwww on 8 Dec 2013
Hi, I am not very good with loops, so I did it be hand for 4 steps and then put them into code, I need somehow to put rFinal into a for loop so goes through 1000 times. I tried to do that, but without success, the loop are very incorrect, but I'm not sure how to fix them Use for reference, the variable status in order: Initial, Next, changed, New, Final Also NOTE: Pi in this case is the velocity not the variable "pi" (had to use variable pi..)
NOTE: the part I need help with is at the bottom (*** after the space), right BEFORE THE FOR LOOP, the rest is just variables if you want to run it.... if true %% Defining the variables
% the masses
M1=606; % galactic mass unit, = 2.32e7 of solar masses % central mass constants
M2=3690; % disk constants
M3=4615; % halo constants
rInitial=8.499; % kpc
rMin=3.51; % kpc
rMax=9.15; % kpc
a2=5.3178; % disk constants % scalar lengths
a3=12; % halo constants % scalar lengths
b1=0.3873; % central mass constants % scalar lengths
b2=0.25; % disk constants % scalar lengths
z=0;
% other remining variables and constants
PiInitial=6.827; % kpc/gtm
dt=1e-3; %time step
%%the order of solving the variables - done by hand
DeltarInitial= PiInitial*dt
rNext= rInitial + DeltarInitial
% solving for rChanged:
% partials of Fi1 % the acceleration components - with r initial
DiyFi1OverRInitial= (M1*rInitial) *dt / ( (b1^2+rInitial^2+z^2)^(3/2) );
% partials of Fi2 % the acceleration components - with r initial
DiyFi2OverRInitial= (M2*rInitial) *dt / ( (rInitial^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
% partials of Fi3 % the acceleration components - with r initial
DiyFi3OverRInitialTerm1= ( 1.02*M3*rInitial*( (sqrt(rInitial^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rInitial^2+z^2)*( 1+( (sqrt(rInitial^2+z^2) / a3)^1.02 )^2 ));
DiyFi3OverRInitialTerm2= ( M3*rInitial*( ( sqrt(rInitial^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rInitial^2+z^2)^(3/2) )*( 1+( ( sqrt(rInitial^2+z^2) / a3 )^1.02 ) ) );
DiyFi3OverRInitialTerm3= ( -2.02*M3*rInitial*( ( sqrt(rInitial^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rInitial^2+z^2)*( 1+( ( sqrt(rInitial^2+z^2) / a3 )^1.02 ) ) ); % a negative term
DiyFi3OverRInitial= DiyFi3OverRInitialTerm1 + DiyFi3OverRInitialTerm2 + DiyFi3OverRInitialTerm3;
% the acceraration component in the r direction
arInitial= DiyFi1OverRInitial + DiyFi2OverRInitial + DiyFi3OverRInitial
PiNext=PiInitial + arInitial*dt
DeltarNext= PiNext*dt
rChanged= rNext + DeltarNext
% sloving for rNew:
% partials of Fi1 % the acceleration components - with rNext
DiyFi1OverRNext= (M1*rNext) *dt / ( (b1^2+rNext^2+z^2)^(3/2) );
% partials of Fi2 % the acceleration components - with rNext
DiyFi2OverRNext= (M2*rNext) *dt / ( (rNext^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
% partials of Fi3 % the acceleration components - with rNext
DiyFi3OverRNextTerm1= ( 1.02*M3*rNext*( (sqrt(rNext^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rNext^2+z^2)*( 1+( (sqrt(rNext^2+z^2) / a3)^1.02 )^2 ));
DiyFi3OverRNextTerm2= ( M3*rNext*( ( sqrt(rNext^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rNext^2+z^2)^(3/2) )*( 1+( ( sqrt(rNext^2+z^2) / a3 )^1.02 ) ) );
DiyFi3OverRNextTerm3= ( -2.02*M3*rNext*( ( sqrt(rNext^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rNext^2+z^2)*( 1+( ( sqrt(rNext^2+z^2) / a3 )^1.02 ) ) ); % a negative term
DiyFi3OverRNext= DiyFi3OverRNextTerm1 + DiyFi3OverRNextTerm2 + DiyFi3OverRNextTerm3;
% the acceraration component in the r direction
arNext= DiyFi1OverRNext + DiyFi2OverRNext + DiyFi3OverRNext
PiChanged=PiNext + arNext*dt
DeltarChanged= PiChanged*dt
rNew= rChanged + DeltarChanged
%
%
%
% sloving for rFinal: - THE LOOP PART******************************************
% partials of Fi1 % the acceleration components - with rChanged
DiyFi1OverRChanged= (M1*rChanged) *dt / ( (b1^2+rChanged^2+z^2)^(3/2) );
% partials of Fi2 % the acceleration components - with rNext
DiyFi2OverRChanged= (M2*rChanged) *dt / ( (rChanged^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
% partials of Fi3 % the acceleration components - with rNext
DiyFi3OverRChangedTerm1= ( 1.02*M3*rChanged*( (sqrt(rChanged^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rChanged^2+z^2)*( 1+( (sqrt(rChanged^2+z^2) / a3)^1.02 )^2 ));
DiyFi3OverRChangedTerm2= ( M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rChanged^2+z^2)^(3/2) )*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) );
DiyFi3OverRChangedTerm3= ( -2.02*M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rChanged^2+z^2)*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) ); % a negative term
DiyFi3OverRChanged= DiyFi3OverRChangedTerm1 + DiyFi3OverRChangedTerm2 + DiyFi3OverRChangedTerm3;
% the acceraration component in the r direction
arChanged= DiyFi1OverRChanged + DiyFi2OverRChanged + DiyFi3OverRChanged
PiNew=PiChanged + arChanged*dt
DeltarNew= PiNew*dt
rFinal= rNew + DeltarNew
for i=1:5
arChanged= DiyFi1OverRChanged + DiyFi2OverRChanged + DiyFi3OverRChanged
arChanged= arChanged + 1
for i=1:5
PiNew(i)=PiChanged + arChanged*dt
for k=1:5
DeltarNew(k)= PiNew*dt
for i=1:5
rFinal(i)= rNew + DeltarNew
end
end
end
end
end
end
Amna on 7 Dec 2013
Sorry fir the late reply: yep the variables I need in order are the ones in the nested loop, but the method (formula) written right above the loop (starting at **) makes more sense, because the i's and k's within the loop make the formula wrong...

sixwwwwww on 7 Dec 2013
The part after *************** can be written in loop as follow:
DiyFi1OverRChanged= (M1 * rChanged) * dt / ( (b1^2+rChanged^2+z^2)^(3/2) );
DiyFi2OverRChanged= (M2*rChanged) *dt / ( (rChanged^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
DiyFi3OverRChangedTerm1= ( 1.02*M3*rChanged*( (sqrt(rChanged^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rChanged^2+z^2)*( 1+( (sqrt(rChanged^2+z^2) / a3)^1.02 )^2 ));
DiyFi3OverRChangedTerm2= ( M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rChanged^2+z^2)^(3/2) )*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) );
DiyFi3OverRChangedTerm3= ( -2.02*M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rChanged^2+z^2)*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) ); % a negative term
DiyFi3OverRChanged= DiyFi3OverRChangedTerm1 + DiyFi3OverRChangedTerm2 + DiyFi3OverRChangedTerm3;
arChanged= DiyFi1OverRChanged + DiyFi2OverRChanged + DiyFi3OverRChanged;
PiNew(1) = PiChanged + arChanged * dt;
DeltarNew(1) = PiNew(1) * dt;
rFinal(1) = rNew + DeltarNew(1);
for i=1:5
arChanged= arChanged + 1;
PiNew(i + 1)= PiChanged + arChanged * dt;
DeltarNew(i + 1)= PiNew(i + 1) * dt;
rFinal(i + 1)= rNew + DeltarNew(i + 1);
end
It it what you need?
sixwwwwww on 8 Dec 2013
Edited: sixwwwwww on 8 Dec 2013
you are welcome.
You can check the values on which rFinal depends. if those values are changing then rFinal must also change accordingly