sciencetalk

Code: Selecteer alles

function gravity

	%gravity at pole=(0,a,0)

	function res=fpole(x,y,z,a)

		res=1/(x^2+(a-y)^2+z^2)*(a-y)/sqrt(x^2+(a-y)^2+z^2);   %|F|*[(pole-r)*1y]/|pole-r|

		if isnan(res)||res<0				 %remove errors due to y~=a

			res=0;

		end

	end



	%gravity at equator (1,0,0)

	function res=feq(x,y,z)

		res=1/((1-x)^2+y^2+z^2)*(1-x)/sqrt((1-x)^2+y^2+z^2); %|F|*[(equator-r)*1x]/|equator-r|

		if isnan(res)||res<0				  %remove errors due to x~=1

			res=0;

		end

	end



grid=0.05;		  %grid for earth

agrid=grid;		 %grid for a

aVector=0:agrid:1;  %vector with values for a

count=0;			%index for vectors

gpole=zeros(1,length(aVector));  %vector for gravity at pole

geq=zeros(1,length(aVector));	%vector for gravity at equator

% hold on;



for a=aVector

	count=count+1;

	for x=-1:grid:1

		for y=-a*sqrt(1-x^2):grid:a*sqrt(1-x^2)   %ellips: x^2+y^2/a^2=1

			for z=-real(sqrt(1-x^2-y^2/a^2)):grid:real(sqrt(1-x^2-y^2/a^2)) %ellipsoid: x^2+y^2/a^2+z^2=1

				gpole(count)=gpole(count)+fpole(x,y,z,a);

				geq(count)=geq(count)+feq(x,y,z);

				%plot3(x,z,y);

			end

		end

	end

end

plot(aVector,gpole,aVector,geq);

end