I was reading about magnetism and it was getting really hard to understand without drawing it out myself. Unfortunately this isn't about magnetism yet, but it IS about potential energy!
Basically, a potential is a field of conservative forces - so if a mass is thrown while subjected to a potential, the value of it's kinetic energy plus its potential energy will remain the same (eg mechanical energy is conserved).
Two great examples of conservative forces are gravity and electrostatics.
The force of gravity is G*m*M*r^(-2). To get the potential energy at some distance r, we integrate with respect to r, giving the potential energy, U=-G*m*M*r^(-1). The potential assumes the object being attracted has a mass of one, so we graph u=-G*M/r. That curve is called the potential. The potential of an object on earth's surface is given by u=g*h. The potential for electrostatics is in the form u=-G*M/r, but it allows for M to be negative.
Plotting this is easy. You have your potential function which is a function of x and y. From x and y (and the position of the mass) you can determine r and plot the function. If you add two masses (so your equation is u=-G*M1/r1-G*M2/r2) you get something like this:
(I take the sine/cosine of u to get pretty colors)
There is however another useful way to visualize potentials: using field lines and equipotentials.
Equipotentials are easy, they're all values of x and y such that u is constant. For the above image, the equipotentials look like this:
Field lines are the opposite. They're the lines going along the path of MAXIMUM change in potential energy. They're also perpendicular to the equipotentials. They look like this:
aand both together:
For fun I reversed the value of M on the larger mass to see what it would look like:
hmm... needs some adjustments:
...Okay, so my program isn't a great artist, but the tech is cool!
Anyways, here's the source code:
global centerx as float
global centery as float
global ppu as float
global w as integer
w=screen width()
global h as integer
h=screen height()
centerx=0
centery=0
ppu=140
for x=0 to w
for y=0 to h
ink colorfunc(potentialfunction(xPixelToUnit(x),yPixelToUnit(y))),0
dot x,y
remstart xu#=xPixelToUnit(x)
yu#=yPixelToUnit(y)
xdu#=xPartialDerivative(xu#,yu#)
ydu#=yPartialDerivative(xu#,yu#)
d#=sqrt(xdu#*xdu#+ydu#*ydu#)
xdu#=xdu#/d#
ydu#=ydu#/d#
line x,y,x-xdu#*5,y+ydu#*5 remend
next y
next x
lock pixels
for n=0 to 10
drawEquipotential(xPixelToUnit(n*w/10.0),yPixelToUnit(h/2))
next n
for n=0 to 15
drawFieldLine(xPixelToUnit(n*w/15.0),yPixelToUnit(0))
next n
for n=0 to 15
drawFieldLine(xPixelToUnit(n*w/15.0),yPixelToUnit(h))
next n
for n=0 to 5
drawFieldLine(xPixelToUnit(0),yPixelToUnit(n*h/5.0))
next n
for n=0 to 5
drawFieldLine(xPixelToUnit(w),yPixelToUnit(n*h/5.0))
next n
unlock pixels
wait key
end
function drawFieldLine(startx as float, starty as float)
n=0 //curent dot number
maxn=20000 //maximum number of integration points
curx#=startx //variables for tracking current pos
cury#=starty
col=rgb(255,255,255)
remstart
singularity checker.
the idea here is that if the integration point has barely moved
in <framestocount> steps, then the calculation will stop. This
happens when points reach a singularity.
remend
framestocount as integer
framestocount=20
dim fx(framestocount-1) as float //0 to framestocount-1
dim fy(framestocount-1) as float
index as integer //current point.
index=1
fx(0)=startx //initialize the array
fy(0)=starty
stepsize#=.001
while n<maxn
//take the normalized gradient
dx#=xPartialDerivative(curx#,cury#)
dy#=yPartialDerivative(curx#,cury#)
d#=sqrt(dx#*dx#+dy#*dy#)
//normalize (dx,dy) and scale by stepsize#.
d#=stepsize#/d#
dx#=dx#*d#
dy#=dy#*d#
//the gradient goes up the potential but we want to go down the potential,
//so subtract (dx,dy) instead of adding it to the current position.
curx#=curx#-dx#
cury#=cury#-dy#
//store it in the singularity checker.
fx(index)=curx#
fy(index)=cury#
inc index //next index
//have index go from 0 to framestocount-1.
if index>=framestocount then index=0
//inc number of points used in integration.
inc n
//draw dots every tenth integration point.
if (n mod 10)=0 then dot xUnitToPixel(curx#),yUnitToPixel(cury#),col
//singularity checker - only starts once the fx/fy arrays have been filled.
if n>=framestocount
subindex=(n+1) mod framestocount //n+1 mod is the OLDEST stored frame in the array.
dx#=curx#-fx(subindex)
dy#=cury#-fy(subindex)
d#=dx#*dx#+dy#*dy# //magnitude difference from the oldest position
if d#<stepsize#*stepsize#*6*6 then n=maxn+1 //if its less than 6*stepsize, exit.
endif
endwhile
endfunction
function drawEquipotential(startx as float, starty as float)
fDrawEquipotential(startx,starty,0)
endfunction
function fDrawEquipotential(startx as float, starty as float, rot as integer)
n=0 //curent dot number
maxn=25000 //maximum number of integration points
stepsize#=.001 //distance from one iteration to another.
//store the current potential for later correction functions.
potential#=potentialFunction(startx,starty)
curx#=startx
cury#=starty
col=rgb(255,255,255)
remstart
with equipotentials you have to check to see if they
loop back on themselves. If they do, the difference between
curx and startx will be near zero. checkstartnum tells the
program to start checking for looping back after checkstartnum
iterations. It can be pretty high, since most equipotential loops are large.
remend
checkstartnum=1000
closedloop=0
while n<maxn
//calculate the gradient vector.
xGradient#=xPartialDerivative(curx#,cury#)
yGradient#=yPartialDerivative(curx#,cury#)
//flip x/y and negate one to rotate 90 degrees.
if rot
dx#=-yGradient#
dy#=xGradient#
else
dx#=yGradient#
dy#=-xGradient#
endif
d#=sqrt(dx#*dx#+dy#*dy#)
d#=d#/stepsize#
dx#=dx#/d#
dy#=dy#/d#
curx#=curx#+dx#
cury#=cury#+dy#
if (n mod 10)=0 then dot xUnitToPixel(curx#),yUnitToPixel(cury#),col
xGradient#=xGradient#/d#
yGradient#=yGradient#/d#
//the difference in potential SHOULD be zero, but due to numerical estimates
//it's not. However, since (yPartial,xPartial) points in the direction of
//increasing potential, we can correct the function easily.
potentialdiff#=potentialFunction(curx#,cury#)-potential#
if potentialdiff#<0
//then we need to increase the current potential:
curx#=curx#+xGradient#
cury#=cury#+yGradient#
else
//else decrease the current potential
curx#=curx#-xGradient#
cury#=cury#-yGradient#
endif
if n>checkstartnum
dx#=curx#-startx
dy#=cury#-starty
d#=dx#*dx#+dy#*dy#
if d#<(10*stepsize#)*(10*stepsize#)
closedloop=1
n=maxn+1
endif
endif
inc n
endwhile
//if its not a closed loop, draw it again in the opposite direction.
if closedloop=0 and rot=0 then fDrawEquipotential(startx,starty,1)
endfunction
function scaleUnitToPixel(u as float)
r#=u
endfunction r#
function scalePixelToUnit(p as float)
r#=p/ppu
endfunction r#
function xUnitToPixel(x as float)
s#=(x-centerx)*ppu+w/2
endfunction s#
function yUnitToPixel(y as float)
s#=-(y-centerx)*ppu+h/2
endfunction s#
//inverses of the above two functions.
function xPixelToUnit(x as float)
s#=(x-w/2)/ppu+centerx
endfunction s#
function yPixelToUnit(y as float)
s#=-(y-h/2)/ppu+centerx
endfunction s#
//randomy gradient color function, for nice colors.
function colorfunc(mag as float)
mag=log(abs(mag)+1) //logarithmic so that inverse square things still look interesting far away.
r#=unitcos(mag*360/3)*255
g#=unitcos(mag*2*360+152)*.1*255
b#=unitcos(mag*2*360)*255
col as dword
col=rgb(int(r#),int(g#),int(b#))
endfunction col
//useful for color functions
function unitcos(n#)
ret#=(cos(n#)+1)/2
endfunction ret#
//takes (pf(x+dx,y)-pf(x,y))/dx where pf is the potentialFunction
function xPartialDerivative(x as float,y as float)
dx#=.01
slope#=(potentialFunction(x+dx#,y)-potentialFunction(x,y))/dx# //numerically differentiate
endfunction slope#
//takes (pf(x,y+dy)-pf(x,y))/dy where pf is the potentialFunction
function yPartialDerivative(x as float,y as float)
dy#=.01
slope#=(potentialFunction(x,y+dy#)-potentialFunction(x,y))/dy# //numerically differentiate
endfunction slope#
//A scalar function of x and y defines a potential.
//CHANGE THIS FUNCTION!
function potentialFunction(x as float, y as float)
d1sq#=(x-1)*(x-1)+y*y
d2sq#=(x+1)*(x+1)+y*y
mag#=-20/sqrt(d1sq#)-10/sqrt(d2sq#)
endfunction mag#
For those of you interested in how it works, you need a bit of calculus. The key term here is GRADIENT. The gradient of u consists of the vector of both partial derivatives of u. so the gradient=(du/dx,du/dy) (partial derivatives). It turns out that this vector is in the direction of the maximum rise in potential. From that you can find the force of the potential (-gradient) and the perpendicular values (gradient rotated by 90 degrees).
