Acabo de re-implementar un viejo programa del Pascal de la Mac del 95.
Adjunto el original y la version para que corra en 3.10.

Hagan en el Workspace 13 timesRepeat: [Planet new].

Calcula por incrementos finitos y alguna vez lo discutimos con Diego Gomez
Deck.

Pensar que cuando lo empece , en la UNRC , le pedi las fórmulas a mis amigos
Físicos y ellos decián (por el 85) que en el observatorio en Córdoba no
tenian para calcular trayectorias de cuerpos celestes para mas de 9 cuerpos.

El artículo que me inspiro lo lei en la version española de Scientific
American.

Tambien en Scientific American me enteré de Xerox Parc, Smalltalk y la
primera versión que salio fuera de Xerox en un contrato de investigación con
Apple, Tektronix, Hewlett Packard y (creo) Digital.

Seguro que puede optimizarse.
Que pasaría si existiese el eter y los planetas tuvieran fricción ?


Edgar

CircleMorph subclass: #Planet instanceVariableNames: 'x y vx vy ax ay m dx dy
cfx cfy mass' classVariableNames: '' poolDictionaries: '' category:
'NewTruch'! !Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! ax
"Answer the value of ax" ^ ax! ! !Planet methodsFor: 'accessing' stamp: 'edc
12/3/2007 11:36'! ax: anObject "Set the value of ax" ax := anObject! !
!Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! ay "Answer the
value of ay" ^ ay! ! !Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007
11:36'! ay: anObject "Set the value of ay" ay := anObject! ! !Planet
methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! cfx "Answer the value of
cfx" ^ cfx! ! !Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'!
cfx: anObject "Set the value of cfx" cfx := anObject! ! !Planet methodsFor:
'accessing' stamp: 'edc 12/3/2007 11:36'! cfy "Answer the value of cfy" ^
cfy! ! !Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! cfy:
anObject "Set the value of cfy" cfy := anObject! ! !Planet methodsFor:
'accessing' stamp: 'edc 12/3/2007 11:36'! dx "Answer the value of dx" ^ dx! !
!Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! dx: anObject "Set
the value of dx" dx := anObject! ! !Planet methodsFor: 'accessing' stamp:
'edc 12/3/2007 11:36'! dy "Answer the value of dy" ^ dy! ! !Planet
methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! dy: anObject "Set the
value of dy" dy := anObject! ! !Planet methodsFor: 'accessing' stamp: 'edc
12/3/2007 11:36'! m "Answer the value of m" ^ m! ! !Planet methodsFor:
'accessing' stamp: 'edc 12/3/2007 11:36'! m: anObject "Set the value of m" m
:= anObject! ! !Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! vx
"Answer the value of vx" ^ vx! ! !Planet methodsFor: 'accessing' stamp: 'edc
12/3/2007 11:36'! vx: anObject "Set the value of vx" vx := anObject! !
!Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007 11:36'! vy "Answer the
value of vy" ^ vy! ! !Planet methodsFor: 'accessing' stamp: 'edc 12/3/2007
11:36'! vy: anObject "Set the value of vy" vy := anObject! ! !Planet
methodsFor: 'stepping and presenter' stamp: 'edc 12/3/2007 19:17'! step | all r
dt p | dt := 1/ self stepTime. all := World submorphsSatisfying: [:morph |
(morph isKindOf: Planet) and: [morph ~= self]]. ax:= ay := dx := dy := 0. x :=
self center x. y := self center y. all do: [:morph | dx := x - morph center x.
dy := y - morph center y. dx := dx * dx. dy := dy* dy. r := (dx + dy)
sqrt. r := r * r . r = 0.0 ifFalse:[ ax := (1000 * mass negated * (x -
morph center x) / r) + ax. ay := (1000 * mass negated * (y - morph center y)
/ r) + ay] ]. vx := vx + (ax * dt). vy := vy + (ay * dt). x := x + (vx * 2 *
dt). y := y + (vy * 2 * dt). p := Point x: x y: y. self center: p. self
comeToFront! ! !Planet methodsFor: 'stepping and presenter' stamp: 'edc
12/3/2007 19:12'! stepTime ^200! ! !Planet methodsFor: 'initialization' stamp:
'edc 12/3/2007 18:46'! initialize | origin extent maxVelocity | super
initialize. maxVelocity := 10. origin _ 0@0. extent _ (World extent - self
bounds extent) rounded. self position: origin x + extent x atRandom @ (origin y
+ extent y atRandom). vx _ 5 atRandom. (2 atRandom = 1) ifTrue:[vx := vx
negated ]. vy _ 5 atRandom. (2 atRandom = 1) ifTrue:[vy:= vy negated ]. ax :=
ay := 0. self color: Color random. mass := 50 atRandom. self extent: self
extent * 5 atRandom. self openInWorld.! !

program Ley; const cuerpos = 13; WindowID = 129; {Resource ID for my
window} var x, y, vx, vy, ax, ay, m, dx, dy: array[1..cuerpos] of real; r,
dt: real; i, j, tiempo, der, aba: integer; myevent: eventrecord; mycolor:
array[1..cuerpos] of RGBColor; ret: rect; screenPort: GrafPtr; myWindow:
WindowPtr; {our window} wRecord: WindowRecord; {storage for
window record} screenRect: rect; {size of screen; could be
machine-dependent} fin: boolean; procedure color_Cuerpos; var i: integer;
begin for i := 2 to cuerpos do begin with mycolor[i] do begin blue
:= random; red := random; green := random; end; end; end; begin
getdatetime(randseed); {Get the port which is the whole screen, to use
when deactivating our window.} { This prevents the current grafPort
pointer from ever dangling.} GetWMgrPort(screenPort); {get whole screen port
that window mgr uses} {get window: use wRecord storage. Port is set to
that of the new window.} { GetNewCWindow posts an update event for
the new window,} { so it will be redrawn right away.} {
This is the one we want on top to start.} myWindow := GetNewCWindow(windowID,
nil, POINTER(-1)); {-1 => frontmost window} SetPort(myWindow); {and
start off with it} ShowWindow(MyWindow); {set up dragRect; we can drag
the window within it} screenRect := screenBits.bounds; {don't assume
screen size} {set drag rect to avoid menu bar, and keep at least 4
pixels on screen} SetRect(ret, 0, 0, screenRect.right, screenRect.bottom); der
:= screenRect.right; aba := screenRect.bottom; dt := 0.5; color_cuerpos; for
i := 2 to cuerpos do begin m[i] := abs(random mod aba) / 5 + 50; x[i] :=
abs(random mod der); y[i] := abs(random mod aba); end; m[1] := der; x[1] :=
der / 2; y[1] := aba / 2; repeat tiempo := tiempo + 1; if tiempo = 1000
then begin eraserect(ret); color_cuerpos; tiempo := 0 end; for
i := 1 to cuerpos do begin rgbForeColor(MyColor[i]);
paintcircle(round(x[i]), round(y[i]), round(m[i] / 10));
invertcircle(round(x[i]), round(y[i]), round(m[i] / 9)); ax[i] := 0; ay[i]
:= 0; for j := 1 to cuerpos do begin if j <> i then begin
dx[i] := (x[i] - x[j]); dx[i] := dx[i] * dx[i]; dy[i] := (y[i] -
y[j]); dy[i] := dy[i] * dy[i]; r := sqrt(dx[i] + dy[i]); if r >
M[i] then begin r := r * r * r; ax[i] := -M[j] * (x[i] - x[j])
/ r + ax[i]; ay[i] := -M[j] * (y[i] - y[j]) / r + ay[i]; end end
end; vx[i] := vx[i] + ax[i] * dt; vy[i] := vy[i] + ay[i] * dt; x[i] :=
x[i] + vx[i] * 2 * dt; y[i] := y[i] + vy[i] * 2 * dt; end; fin :=
GetNextEvent(everyEvent, myEvent); until (myEvent.what = keydown) or
(myEvent.what = mousedown) end.