\ ob.gravity
\
\ hp_gravity
\ class of force that acts like the thing that dropped the apple.
\ part of the collection of software componenets: "henri poincare."
\
\ based on ideas presented by Joel Ryan.
\
\ code: Han-earl Park
\ copyright 2004 buster & friends' C-ALTO Labs
\ (Den haag, December 1997 -
\ (Valencia, October 1998 -
\ (Southampton, February 2002 -
\
\ MOD: HeP 02/11/98	Started gravity class.
\ MOD: HeP 11/30/98	Change the stack diagram of the (calc.force): method
\					from (  -- fx fy fz... ) to (  --  ) which should make
\					the method easier to code.
\ MOD: HeP 07/27/99	Change name of components to "henri poincare"
\ MOD: HeP 10/26/99	Finally add ability for gravitational force to act
\					across toroid shape. The current implementation is
\					a mess, but it works!
\ MOD: HeP 02-23-02	Use the force constant from the OB.FORCE class. No longer
\					have the specialized iv-grav-constant.
\ MOD: HeP 03-19-04	Fixed the kludgy stack behavior in ABS<.
\ MOD: HeP 03-20-04	Cleanup code GRAV.CALC.DISP calls GRAV.CALC.DIM.DISP in
\					its loop.
\					Remove redundant calcualtions when parallel updating of
\					particles.
\
\ ToDo: store vectored displacement between particles internally
\ ToDo:	NO internal scaling...?
\		Use double numbers....?
\		Use floating point....?

include?  task-hp_force  myt:hp_force


anew task-hp_gravity


:class  OB.GRAVITY  <super  OB.FORCE
	
	
:m DEFAULT:	(  --  )
	default: super
	40 iv=> iv-force-constant
;m


: ABS<	( n1 n2 -- flag )
	ABS  swap  ABS
	>
;

: ABS.MIN	( n1 n2 -- min , absolute min )
	2dup  ABS<
	IF	drop  ELSE  nip
	THEN
;

: GRAV.CALC.DIM.DISP	{ p1 p2 dim# -- disp }
	hp_position dim#  ed.at: p2
	hp_position dim#  ed.at: p1  -
	\
	dim# get.wrap: p1
	IF
		hp_position dim#  ed.at: p2	\ x2
		hp_position dim#  ed.at: p1	\ x1
		-									\ x2 - x1
		dim# get.range: p1  NIP				\ range
	\
		2dup  -								\ x2 -  x1 - range
		-rot  +								\ x2 -  x1 + range
	\
		ABS.MIN
	\
		ABS.MIN
	THEN
;

: GRAV.CALC.DISP	( p1 p2 -- disp , calculate displacement )
	#dimensions@ 0
	DO
		2dup i	GRAV.CALC.DIM.DISP
		-rot
	LOOP
	2drop								
	vec->abs
;

: GRAV.COLLISION?	{ p1 p2 disp -- flag , coliision of particles? }
	get.radius: p1
	get.radius: p2
	+
	disp abs
	>
;

: GRAV.CALC.FORCE	( disp m1 m2 -- force , gravitational pull )
	*					(  -- disp m1*m2 )
	iv-force-constant *	(  -- disp m1*m2*c )
	force.scale.value	(  -- disp m1*m2*c' )
	swap				(  -- m1*m2*c disp )
\
	#dimensions@ 2
	DO
		dup  *
	LOOP
	/
;

\ force vectors ar proportional to spacial displacement:
\
\	fx = F * x / disp
\

: GRAV.UPDATE.FORCE	{ p1 p2 disp force --  }
\	force  force.unscale.value  -> force
\
	#dimensions@ 0
	DO
		p1 p2 i GRAV.CALC.DIM.DISP
	\
		force *
		disp /mod		(  -- mod fx )
	\
		dup
		IF
			nip					(  -- fx )
			force.unscale.value	(  -- fx )
		ELSE
			drop
			IF  1  ELSE -1
			THEN
		THEN
	\
		hp_force i  ed.at: p1  +	\ append value to current
		hp_force i  ed.to: p1
	LOOP
;

:m (ELM.CALC.FORCE):	(  --  )
	iv-force-#pucks 0
	DO
		i iv-force-cur-puck#  =  not	\ ignore if current particle
		IF
			iv-force-cur-puck				(  -- p1 )
			i iv-force-cur-space  get: []	(  -- p1 p2 )
		\
			2dup  GRAV.CALC.DISP			(  -- p1 p2 disp )
			3dup							(  -- p1 p2 disp p1 p2 disp )
		\
			-rot  get.mass: []				(  -- p1 p2 disp disp p1 m2 )
			swap  get.mass: []				(  -- p1 p2 disp disp m2 m1 )
			GRAV.CALC.FORCE					(  -- p1 p2 disp force )
		\
			GRAV.UPDATE.FORCE				(  --  )
		THEN
	LOOP
;m


\ update particles in parallel

: PARALLEL.GRAV.UPDATE.FORCE	{ p1 p2 disp force --   }
	#dimensions@ 0
	DO
		p1 p2 i GRAV.CALC.DIM.DISP		\  -- dim_disp
	\
		force *							\  -- dim_disp*force
		disp /							\  -- dim_force
	\
		dup
		hp_force i  ed.at: p1  +		\ sum value to current
		hp_force i  ed.to: p1			\ update
	\
		negate
		hp_force i  ed.at: p2  +		\ sum value to current
		hp_force i  ed.to: p2			\ update
	LOOP
;

: PARALLEL.GRAV.CALC.FORCE	(  --  )
	iv-force-#pucks iv-force-cur-puck# 1+	\ update only those needed
	DO
		iv-force-cur-puck				\  -- p1
		i iv-force-cur-space  get: []	\  -- p1 p2
	\
		2dup GRAV.CALC.DISP				\  -- p1 p2 disp
		3dup							\  -- p1 p2 disp p1 p2 disp
	\
		swap  get.mass: []				\  -- p1 p2 disp p1 disp m2
		rot   get.mass: []				\  -- p1 p2 disp disp m2 m1
		GRAV.CALC.FORCE					\  -- p1 p2 disp force )
	\
		PARALLEL.GRAV.UPDATE.FORCE		\  --
	LOOP
;

:m (CALC.FORCE):	(  --  , calculate forces for all particles )
	iv-force-#pucks 0
	DO
		i iv-force-cur-space get: []			\  -- p
	\
		dup  get.type: []  self get.type: []	\  -- p t1 t2
		particle.type=
		IF
			   iv=> iv-force-cur-puck
			i  iv=> iv-force-cur-puck#
		\
			PARALLEL.GRAV.CALC.FORCE			\  --
		ELSE
			drop
		THEN
	LOOP
;m

;class