; physlaws.lsp              Gordon S. Novak Jr.            ; 25 Aug 06

; Copyright (c) 2006 Gordon S. Novak Jr. and The University of Texas at Austin.

; Definitions for physical laws.  Derived from equations.lsp

; 19 Feb 04; 23 Feb 04; 27 Feb 04; 01 Mar 04; 02 Mar 04; 15 Mar 04; 17 Mar 04
; 18 Mar 04; 19 Mar 04; 22 Mar 04; 23 Mar 04; 30 Mar 04; 31 Mar 04; 01 Apr 04
; 02 Apr 04; 05 Apr 04; 06 Apr 04; 08 Apr 04; 09 Apr 04; 12 Apr 04; 13 Apr 04
; 15 Apr 04; 16 Apr 04; 20 Apr 04; 22 Apr 04; 23 Apr 04; 27 Apr 04; 03 Jun 04
; 07 Apr 05; 13 Sep 05; 06 Dec 05; 28 Jun 06; 29 Jun 06; 30 Jun 06

; Things that should be added:
;   symbol used externally for each variable (e..g. rho for resistivity)

(defmacro basis-vars     (x) `(get ,x 'basis-vars))
(defmacro vars-units     (x) `(get ,x 'vars-units))
(defmacro vars-unittypes (x) `(get ,x 'vars-unittypes))
(defmacro equations      (x) `(get ,x 'equations))
(defmacro assumed-unit   (x) `(get ,x 'assumed-unit))

(defvar *assumed-units*)
(setq *assumed-units* '(
  ((time period) second)
  (frequency hertz)
  ((length radius height distance wavelength) meter)
  (mass kilogram)
  ((velocity speed) (/ meter second))
  ((acceleration gravity) (/ meter (* second second)))
  (area (* meter meter))
  (volume (* meter meter meter))
  ((weight force) newton)
  (pressure pascal)
  ((impulse momentum) (* kilogram (/ meter second)))
  ((energy kinetic-energy potential-energy work) joule)
  (power watt)
  (charge coulomb)
  (capacitance farad)
  ((voltage emf) volt)
  (magnetic-field tesla)
  (resistance ohm)
  (current ampere)
  (cost dollar)
  (temperature degree-kelvin)
  ))

(dolist (pair *assumed-units*)
  (if (consp (car pair))
      (dolist (var (car pair)) (setf (assumed-unit var) (cadr pair)))
      (setf (assumed-unit (car pair)) (cadr pair)) ) )

; format is (name basis-vars constant-vars units equations)
; units is (var unit assumed-value) where assumed-value is used in changeobjvar
(defvar *physlaws*)
(setq *physlaws* '(
; Geometry

(circle (radius centerx centery)
  ((diameter meter)
   (circumference meter))
  ((= diameter (* 2 radius))
   (= circumference (* pi diameter))
   (= area (* pi (expt radius 2))) ) )

(sphere (radius)
  ((diameter meter)
   (circumference meter))
  ((= diameter (* 2 radius))
   (= circumference (* pi diameter))
   (= area (* 4 (* pi (expt radius 2))))
   (= volume (* (/ (* 4 pi) 3) (expt radius 3)))
   (= center (tuple (x centerx) (y centery)))
   (= centerx (x center))
   (= centery (y center)) ) )

(cylinder (radius length)
  ((diameter meter)
   (circumference meter))
  ((= diameter (* 2 radius))
   (= circumference (* pi diameter))
   (= area (* circumference length))
   (= volume (* (* pi (expt radius 2)) length)) ) )

(cone (radius height)
  ((side meter)
   (side-area (* meter meter))
   (bottom-area (* meter meter))
   (total-area (* meter meter)) )
  ((= diameter (* 2 radius))
   (= circumference (* pi diameter))
   (= side (sqrt (+ (expt radius 2) (expt height 2))))
   (= side-area (* pi (* radius side)))
   (= bottom-area (* pi (expt radius 2)))
   (= volume (* (/ pi 3) (* (expt radius 2) height)))
   (= total-area (+ side-area bottom-area)) ))

(square (side)
  ((side meter)
   (circumference meter)
   (diagonal meter))
  ((= diagonal (* (sqrt 2) side))
   (= circumference (* 4 side))
   (= area (expt side 2)) ) )

(rectangle (width height)
  ((width meter)
   (height meter)
   (circumference meter)
   (diagonal meter))
  ((= diagonal (sqrt (+ (expt width 2) (expt height 2))))
   (= circumference (* 2 (+ width height)))
   (= area (* width height)) ) )

; density and number of items in area/volume assuming perfect packing
(conglomerate (number item-volume)
  ((number unity))
  ((= number (/ container-volume item-volume))))


; 01 Mar 96
(line-segment (p1x p1y p2x p2y)
  ()
  ((= p1     (tuple (x p1x) (y p1y)))
   (= p1x    (x p1))
   (= p1y    (y p1))
   (= p2     (tuple (x p2x) (y p2y)))
   (= p2x    (x p2))
   (= p2y    (y p2))
   (= deltax (- p2x p1x))
   (= deltay (- p2y p1y))
   (= size   (tuple (x deltax) (y deltay)))
   (= deltax (x size))
   (= deltay (y size))
   (= slope  (/ deltay (float deltax)))
   (= slope  (tan theta))
   (= slope  (/ 1.0 (tan phi)))
   (= length (sqrt (+ (expt deltax 2) (expt deltay 2))))
   (= theta  (atan deltay deltax))
   (= phi    (- (/ pi 2.0) theta))
   (= phi    (atan deltax deltay))
   (= deltay (* length (sin theta)))
   (= deltax (* length (cos theta)))
   (= deltay (* length (cos phi)))
   (= deltax (* length (sin phi))) ) )

(vector (x y)
  ()
  ((= slope  (/ y (float x)))
   (= slope  (tan angle))
   (= slope  (/ 1.0 (tan phi)))
   (= magnitude (sqrt (+ (expt x 2) (expt y 2))))
   (= angle  (atan y x))
   (= phi    (- (/ pi 2.0) angle))
   (= phi    (atan x y))
   (= y      (* magnitude (sin angle)))
   (= x      (* magnitude (cos angle)))
   (= y      (* magnitude (cos phi)))
   (= x      (* magnitude (sin phi))) ) )

(region (left bottom width height)
  ((left meter) (bottom meter) (width meter) (height meter))
  ((= left   (x start))
   (= Bottom (y start))
   (= start  (tuple (x left) (y bottom)))
   (= width  (x size))
   (= height (y size))
   (= area   (* width height))
   (= perimeter (+ (* 2 width) (* 2 height)))
   (= size   (tuple (x width) (y height)))
   (= right  (+ left width))
   (= top    (+ bottom height))
   (= c2     (tuple (x right) (y bottom)))
   (= right  (x c2))
   (= bottom (y c2))
   (= c3     (tuple (x right) (y top)))
   (= right  (x c3))
   (= top    (y c3))
   (= c4     (tuple (x left) (y top)))
   (= left   (x c4))
   (= top    (y c4)) ) )

(right-triangle (x y)
  ((x meter) (y meter))
  ((= theta  (atan y x))
   (= phi    (atan x y))
   (= r      (sqrt (+ (expt x 2) (expt y 2))))
   (= x      (* r (cos theta)))
   (= y      (* r (sin theta)))
   (= x      (* r (sin phi)))
   (= y      (* r (cos phi)))
   (= slope  (tan theta))
   (= slope  (/ y x)) ) )

(triangle (side-a side-b side-c)
  ((side-a meter) (side-b meter) (side-c meter))
  ((= angle-a (- '(q 180 degrees) (+ angle-b angle-c)))
   (= side-a (* side-b (/ (sin angle-a) (sin angle-b))))
   (= side-b (* side-a (/ (sin angle-b) (sin angle-a))))
   (= side-c (* side-b (/ (sin angle-c) (sin angle-b))))
   (= angle-a (acos (/ (- (+ (expt side-b 2) (expt side-c 2))
			  (expt side-a 2))
		       (* 2 (* side-b side-c))))) ) )

; for block on inclined plane problems
(inclined-plane (x y weight)
  ((r meter) (x meter) (y meter) (theta radian) (phi radian)
   (slope unity) (normal-force newton) (tangential-force newton)
   (friction-force newton) (net-force newton) (mu unity) )
  ((= gravity     '(q 9.80665 (/ m (* s s))))
   (= mu 0)   ; default
   (= weight      (* gravity mass))
   (= theta  (atan y x))
   (= phi    (atan x y))
   (= r      (sqrt (+ (expt x 2) (expt y 2))))
   (= x      (* r (cos theta)))
   (= y      (* r (sin theta)))
   (= x      (* r (sin phi)))
   (= y      (* r (cos phi)))
   (= slope  (tan theta))
   (= slope  (/ y x))
   (= normal-force (* weight (cos theta)))
   (= tangential-force (* weight (sin theta)))
   (= friction-force (* mu normal-force))
   (= work    (* friction-force r))
   (= net-force (- tangential-force friction-force))
   (= net-force (* mass acceleration))
   (= r (* 1/2 (* acceleration (expt time 2))))
   (= velocity (* acceleration time))
 ))

; Mechanics

; perhaps should have a 'mechanics' task for e.g. Princeton 71.

(uniform-motion (x0 v dt)
  ((x0 meter) (distance meter) (velocity (/ meter second))
   (dt second) (t0 second))
  ((= x (+ x0 distance))
   (= t0 '(q 0 second))
   (= x0 '(q 0 meter))
   (= dt (- time t0))
   (= distance (* velocity dt)) ))

(acceleration (force mass)
  ()
  ((= force  (* mass acceleration)) ) )   ; Newton's 2nd law

(uniform-acceleration (acceleration time)
  ((distance meter) (acceleration (/ meter (* second second)))
   (velocity (/ meter second)))
  ((= distance      (* (/ acceleration 2) (expt time 2)))
   (= velocity      (* acceleration time)) ) )

(falling (time)
  ((g (/ meter (* second second))) (h meter) (v (/ meter second)))
  ((= g      '(q 9.80665 (/ m (* s s))))
   (= h      (* (/ g 2) (expt time 2)))
   (= v      (* g time)) ) )

(fall (time)
  ((horizontal-distance meter)
   (horizontal-velocity (/ meter second))
   (total-velocity (/ meter second)))
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= horizontal-velocity '(q 0 (/ m s)))           ; default
   (= height   (* 1/2 (* gravity (expt time 2))))
   (= velocity (* gravity time))                   ; vertical velocity
   (= kinetic-energy (* 1/2 (* mass (expt total-velocity 2))))
   (= horizontal-distance (* horizontal-velocity time))
   (= total-velocity (sqrt (+ (expt velocity 2) (expt horizontal-velocity 2))))
  ))

; (phys '(what is the distance of a projectile with height 7.5 m and
;         x-velocity 4.5 m/s))    ; Giancoli 3.19 p. 72
; (phys '(what is the x-velocity of a projectile with height 56 m and
;         distance 45 m))         ; Giancoli 3.24 p. 72
; angle is w.r.t. horizontal
; average acceleration, impulse, time of max height   AP #21
(projectile (time)
  ((angle unity) (x-velocity (/ meter second)) (y-velocity (/ meter second))
   (v0 (/ meter second)))
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= v1 v0)
   (= x-velocity (* v0 (cos angle)))
   (= y-velocity (* v0 (sin angle)))
   (= time (/ (* 2 y-velocity) gravity))     ; time to impact at y = 0
   (= distance (* time x-velocity))          ; x distance to impact at y = 0
   (= height   (/ (expt y-velocity 2) (* 2 gravity)))
 ))

(centrifugal-force (m v r)
  ((acceleration (/ meter (* second second))) (f newton) (v (/ meter second))
   (r meter))
  ((= acceleration (/ (expt v 2) r))
   (= f            (* m acceleration))))

(circular-motion (mass velocity radius)

  ((omega (/ radian second))
   (moment-of-inertia (* kilogram (* meter meter)))
   (angular-momentum (/ (* kilogram (* meter meter)) second)) )
  ((= acceleration (/ (expt velocity 2) radius))
   (= force        (* mass acceleration))
   (= kinetic-energy   (* (/ mass 2) (expt velocity 2)))
   (= moment-of-inertia (* mass (expt radius 2)))
   (= omega (/ velocity radius))
   (= angular-momentum (* omega moment-of-inertia)) ))

(gravity-kinetic (mass height velocity)
  ()
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= kinetic-energy   (* (/ mass 2) (expt velocity 2)))
   (= potential-energy (* (* mass gravity) height))
   (= energy (+ kinetic-energy potential-energy))) )

(kinetic-energy (mass velocity)
  ()
  ((= energy (* (/ mass 2) (expt velocity 2))) ))

(work-energy (mass v1 v2)
  ((v1 (/ meter second)) (v2 (/ meter second))
   (energy1 joule) (energy2 joule))
  ((= energy-change (- energy2 energy1))
   (= work energy-change)
   (= power (/ work time))
   (= energy1 (* 1/2 (* mass (expt v1 2))))
   (= energy2 (* 1/2 (* mass (expt v2 2))))
   ))

; ***** add to this.  AP #22
; add frequency of oscillation of spring Princeton #71
(harmonic-motion (mass height velocity)
  ((angular-speed (/ radian second)))
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= kinetic-energy   (* (/ mass 2) (expt velocity 2)))
   (= potential-energy (* (* mass gravity) height))
   (= energy (+ kinetic-energy potential-energy))
   (= speed (* frequency wavelength))
   (= angular-speed (* (* 2 pi) frequency))
   (= period (/ 1 frequency)) ))

; These equations assume that m1 is the large body (e.g. planet)
; and m2 is the small body (e.g. satellite).
(gravitation (m1 m2 radius)
  ((m1 kilogram) (m2 kilogram) (escape-velocity (/ meter second)))
  ((= gravitation-constant '(q 6.6742e-11 (/ (* nt m m) (* kg kg))))
   (= gravity (/ (* gravitation-constant m1) (expt radius 2)))
   (= satellite-velocity (sqrt (/ (* m1 gravitation-constant) radius)))
   (= kinetic-energy   (* (/ m2 2) (expt satellite-velocity 2))) ; satellite
   (= potential-energy (- (/ (* gravitation-constant (* m1 m2)) radius)))
   (= escape-velocity (sqrt (/ (* 2 (* m1 gravitation-constant)) radius)))
   (= force   (/ (* gravitation-constant (* m1 m2)) (expt radius 2)))
   (= moment-of-inertia (* m2 (expt radius 2)))                  ; satellite
   (= omega (/ satellite-velocity radius))                       ; satellite
   (= angular-momentum (* omega moment-of-inertia))              ; satellite
 ))

(weight (mass volume)
  ()
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= weight  (* gravity mass))
   (= density (/ mass volume))  ))

(physob (mass volume)
  ()
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= weight  (* gravity mass))
   (= density (/ mass volume))  ))

(friction (normal-force mu)
  ((friction-force newton) (mu unity) (normal-force newton))
  ((= friction-force (* mu normal-force))
   (= work    (* friction-force distance))
   (= power   (* friction-force speed)) ))

(friction-weight (weight mu)
  ((friction-force newton) (mu unity))
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= friction-force (* mu weight))
   (= work    (* friction-force distance))
   (= power   (* friction-force speed))
   (= weight  (* gravity mass)) ))

; two masses, one on table (m1), one connected via string and pulley (m2),
; possibly with friction
(weight-pulley (m1 m2)
  ((m1 kilogram) (m2 kilogram) (friction-force newton) (net-force newton)
   (m1-weight newton) (m2-weight newton) (total-mass kilogram) (mu unity))
  ((= mu 0)                     ; default
   (= gravity '(q 9.80665 (/ m (* s s))))
   (= friction-force (* mu m1-weight))
   (= work    (* friction-force distance))
   (= power   (* friction-force speed))
   (= m1-weight  (* gravity m1))
   (= net-force (- m2-weight friction-force))
   (= m2-weight (* gravity m2))
   (= net-force (* total-mass acceleration))
   (= total-mass (+ m1 m2))
 ))

(spring (spring-constant stretch)
  ((spring-constant unity) (stretch meter))
  ((= gravity '(q 9.80665 (/ m (* s s))))
   (= weight  (* gravity mass))
   (= force weight)
   (= force (* spring-constant stretch)) ))   ; must be careful about direction

(lift (mass height time)
  ()
  ((= gravity     '(q 9.80665 (/ m (* s s))))
   (= weight      (* gravity mass))
   (= acceleration '(q 0 (/ m (* s s))))  ; default.  acceleration upward
   (= force (+ weight (* mass acceleration)))
   (= work        (* force height))
   (= speed       (/ height time))
   (= power       (* force speed))
   (= power       (/ work time)) ))

(collide-and-stick (m1 v1 m2 v2)
  ((total-mass kilogram) (total-momentum (* kilogram (/ meter second)))
   (m1 kilogram) (m2 kilogram) (v1 (/ meter second)) (v2 (/ meter second))
   (final-velocity (/ meter second)))
  ((= total-mass (+ m1 m2))
   (= total-momentum (+ (* m1 v1) (* m2 v2)))
   (= total-momentum (* total-mass final-velocity))))

(momentum (mass velocity)
  ()
  ((= momentum (* mass velocity))
   (= impulse momentum)))

; two bodies exert a force on each other.
(symmetric-forces (m1 m2 force)
  ((m1 kilogram) (m2 kilogram) (f1 newton) (f2 newton)
   (a1 (/ meter (* second second))) (a2 (/ meter (* second second))) )
  ((= f1 force)
   (= f2 (- f1))
   (= f1 (* m1 a1))
   (= f2 (* m2 a2))
   ))

(rope (tension)
  ((tension newton))
  ((= force (- tension))))    ; tension >= 0

; object floating in a liquid
(float (mass volume)
  ((displaced-volume (* meter meter meter))
   (floating-volume (* meter meter meter))
   (percentage-floating unity)
   (percentage-below unity)
   (density (/ kilogram (* meter meter meter))))
  ((= gravity     '(q 9.80665 (/ m (* s s))))
   (= weight      (* gravity mass))
   (= density  (/ weight volume))
   (= liquid-density '(q 1000 (/ kilogram (* meter meter meter))))   ; default: water
   (= displaced-volume (/ mass liquid-density))
   (= floating-volume (- volume displaced-volume))
   (= proportion-floating (/ floating-volume volume))
   (= proportion-below (/ displaced-volume volume))
 ))

; determine density of an object by weighing it in fluid and out
; used as an example for entering new knowledge
(object-density (weight mass)
  ((displaced-volume (* meter meter meter))
   (weight-in-liquid newton) (displaced-weight newton)
   (displaced-mass kg) )
  ((= gravity '(q 9.80665 (/ m (* s s))))           ; default
   (= weight  (* gravity mass))
   (= mass    (* density volume))
   (= weight-in-liquid (- weight displaced-weight))
   (= displaced-weight (* gravity displaced-mass))
   (= fluid-density '(q 1 (/ kg liter)))            ; water, default
   (= displaced-mass (* fluid-density volume))
 ))


; Electronics

(electric-dipole (q1 q2 radius)
  ((electric-field (/ volt meter)) (q1 coulomb) (q2 coulomb))
  ((= ke '(q 8.9875517879979115E9
	     (/ (* Newton (* meter meter)) (* coulomb coulomb))))
   (= electric-field (* ke (/ q1 (expt radius 2))))
   (= force (* q2 electric-field))
   (= potential-energy (* ke (* q2 (/ q1 radius))))
   (= work potential-energy)
 ))

(charged-sphere (q radius)
  ((electric-field (/ volt meter)) (q coulomb) (potential volt))
  ((= ke '(q 8.9875517879979115E9
	     (/ (* Newton (* meter meter)) (* coulomb coulomb))))
   (= electric-field (* ke (* q (expt radius 2))))
   (= potential (/ (* ke q) radius))
  ))

(capacitor (capacitance voltage)
  ((electric-field (/ volt meter)))
  ((= epsilon0 '(q 8.854187817E-12 (/ farad meter)))
   (= electric-field (/ voltage distance))
   (= force (* charge electric-field))
   (= dielectric-constant 1)  ; default
   (= capacitance (* (* epsilon0 dielectric-constant)
		     (/ area distance)))
   (= charge (* capacitance voltage))
   (= work energy)
   (= energy (* 1/2 (* capacitance (expt voltage 2))))
  ))

; *****
; add equation for induced current 1998 #45
(magnetic-field (resistance dflux/dt)
  ((delta-time second) (field tesla))
  ((= flux (* field area))
   (= dflux/dt (/ delta-flux delta-time))
   (= emf dflux/dt)
   (= emf (* current resistance))     ; induced current in a wire loop
))


; *****
; This combines (a) emf due to wire moving in magnetic field
;               (b) magnetic field due to current in a wire
;               (c) force on a wire in a magnetic field
; Maybe they should be separate
(wire-magnetic (length magnetic-field)
  ()
  ((= mu0 '(q 12.566370614359172E-7))  ; exactly 4 * PI * 10^-7
   (= emf (* length (* magnetic-field velocity)))
   ; force/length in a current-carrying wire, assumed perpendicular to field
   (= force/length (* magnetic-field current))
   (= magnetic-field (/ (* mu0 current) (* (* 2 pi) radius)))
	))

(parallel-capacitors (c1 c2)
  ((c1 farad) (c2 farad))
  ((= capacitance (+ c1 c2)) ))

(series-capacitors (c1 c2)
  ((c1 farad) (c2 farad))
  ((= capacitance  (/ 1 (+ (/ 1 c1) (/ 1 c2)))) ))

(resistor (resistance current)
  ((unit-cost (/ dollar joule)))
  ((= voltage (* resistance current))
   (= power   (* resistance (expt current 2)))
   (= power   (/ (expt voltage 2) resistance))
   (= work    (* power time))
   (= cost    (* work unit-cost)) ))

; ohmic conductor
(wire (resistivity radius length)
  ((resistivity (* ohm meter)) (diameter meter) (circumference meter)
   (resistance/length (/ ohm meter)))
  ((= diameter (* 2 radius))
   (= circumference (* pi diameter))
   (= area (* pi (expt radius 2)))
   (= resistance/length (/ resistivity area))
   (= resistance (* resistance/length length)) ))

(series-resistors (r1 r2)
  ((r1 ohm) (r2 ohm))
  ((= resistance (+ r1 r2)) ))

(parallel-resistors (r1 r2)
  ((r1 ohm) (r2 ohm))
  ((= resistance (/ 1 (+ (/ 1 r1) (/ 1 r2)))) ))

(battery (voltage current)
  ((loss-voltage volt) (loss-power watt) (terminal-voltage volt)
   (internal-resistance ohm))
  ((= loss-voltage (* internal-resistance current))
   (= loss-power (* internal-resistance (expt current 2)))
   (= emf voltage)
   (= terminal-voltage (- voltage loss-voltage))
   (= power (* terminal-voltage current))
   (= work (* charge terminal-voltage))
   (= internal-resistance '(q 0 ohm))
 ))

(fluid-flow (area velocity)
  ((diameter meter) (circumference meter)
   (c2 meter) (d2 meter) (area2 (* meter meter))
   (flow-rate (/ (* meter meter meter) second)))
  ((= diameter (* 2 radius))
   (= d2 (* 2 r2))
   (= circumference (* pi diameter))
   (= c2 (* pi d2))
   (= area (* pi (expt radius 2)))
   (= area2 (* pi (expt r2 2)))
   (= flow-rate (* velocity area))
   (= flow-rate (* v2 area2))
   (= bogo-bernoulli (* velocity pressure))  ; qualitative only
 ))

; waves

; Assumes wavelength = distance between 2 nodes in a standing wave
(wave (frequency speed)
  ((fundamental-frequency (/ 1 second)) (angular-speed (/ radian second))
   (period second))
  ((= speed (* frequency wavelength))
   (= fundamental-frequency (/ frequency 2))
   (= angular-speed (* (* 2 pi) frequency))
   (= period (/ 1 frequency)) ))
 
; ***** relate mass/length, tension, etc.  1998 #FR5
; harmonic-frequency = frequency of the nth harmonic
(string-wave (fundamental-frequency mass/length length)
  ((frequency hertz) (harmonic-frequency hertz) (harmonic-number unity))
  ((= fundamental-frequency (/ (sqrt (/ force mass/length)) (* 2 length)))
             ; Serway p. 443
   (= harmonic-frequency (* harmonic-number fundamental-frequency))
   ))

(diffraction (wavelength offset-first-max screen-distance)
  ((offset-first-max meter) (screen-distance meter)
   (angle radian) (slit-separation meter))
  ((= offset-first-max (* (sin angle) radius))
   (= radius (sqrt (+ (expt offset-first-max 2)
		      (expt screen-distance 2))))
   (= wavelength (* slit-separation (/ offset-first-max radius)))
   (= speed (* frequency wavelength))
   (= period (/ 1 frequency)) ))

(sound (frequency)
  ((doppler-frequency hertz)
   (observer-speed (/ meter second))
   (source-speed (/ meter second))
   (intensity (/ watt (* meter meter)))
   (power watt))
  ((= speed '(q 343 (/ meter second)))   ; air at 20 C
   (= speed (* frequency wavelength))
   (= period (/ 1 frequency))
   (= doppler-frequency (* frequency (/ (+ speed observer-speed)
					(- speed source-speed))))
   (= intensity0 '(q 1.0e-12 (/ watt (* meter meter))))
   (= decibel (* 10 (log10 (/ intensity intensity0))))
   (= intensity (/ power (* (* 4 pi) (expt radius 2))))
 ))

(light (frequency)
  ((index-of-refraction unity))
  ((= speed-of-light '(q 299792458 (/ meter second)))
   (= planck-constant '(q 6.6260693E-34 (* joule second)))
   (= index-of-refraction 1)     ; default
   (= speed (/ speed-of-light index-of-refraction))
   (= speed (* frequency wavelength))
   (= momentum (/ planck-constant wavelength))
   (= period (/ 1 frequency))
   (= energy (* planck-constant frequency))
 ))

(refraction (index-of-refraction1 index-of-refraction2)
  ((index-of-refraction1 unity) (index-of-refraction2 unity)
   (theta1 radian) (theta2 radian))
  ((= speed-of-light '(q 299792458 (/ meter second)))
   (= speed (/ speed-of-light index-of-refraction))
   (= speed (* frequency wavelength))
   (= period (/ 1 frequency))
	(= (* index-of-refraction1 (sin theta1))
	   (* index-of-refraction2 (sin theta2)))
 ))

(concave-mirror (focal-length)
  ((focal-length meter 0.1)   ; to avoid divide by 0 in changeobjvar
   (image-distance meter)
   (subject-distance meter) (image-height meter) (subject-height meter)
   (magnification unity))
  ((= focal-length (/ radius 2))
   (= (/ 1 focal-length) (+ (/ 1 image-distance) (/ 1 subject-distance)))
   (= magnification (- (/ image-distance subject-distance)))
   (= image-height (* magnification subject-height))
 ))

(converging-lens (focal-length)
  ((focal-length meter) (magnification unity)
   (image-distance meter) (subject-distance meter)
   (image-height meter) (subject-height meter))
  ((= focal-length (/ radius 2))
   (= (/ 1 focal-length) (+ (/ 1 image-distance) (/ 1 subject-distance)))
   (= magnification (- (/ image-distance subject-distance)))
   (= image-height (* magnification subject-height))
 ))

(photoelectric (frequency)
  ((work-function electron-volt))
  ((= speed-of-light '(q 299792458 (/ meter second)))
   (= planck-constant '(q 6.6260693E-34 (* joule second)))
   (= speed-of-light (* frequency wavelength))
   (= work-function '(q 4.08 electron-volt))     ; default: Al
   (= energy (* planck-constant frequency))
   (= kinetic-energy (- energy work-function))
   (= fudge-factor '(q 1 (/ ampere candela)))   ; ***
   (= current (* fudge-factor intensity))       ; *** approx
   (= electron-mass '(q 9.1093826E-31 kg))
   (= kinetic-energy (* 1/2 (* electron-mass (expt velocity 2))))
 ))

; particle physics

(elementary-particle ()
  ()
  ((= planck-constant '(q 6.6260693E-34 (* joule second)))
   (= momentum (* mass velocity))
   (= wavelength (/ planck-constant momentum))
   (= kinetic-energy (* 1/2 (* mass (expt velocity 2))))
 ))

(radioactive-decay (initial-rate half-life)
  ((initial-rate hertz) (final-rate hertz) (half-life second)
   (initial-amount kilogram) (final-amount kilogram))
  ((= final-rate (/ initial-rate (expt 2 (/ time half-life))))
   (= final-amount (/ initial-amount (expt 2 (/ time half-life))))
 ))

(nuclear-reaction (final-mass initial-mass)
  ((final-mass kilogram) (initial-mass kilogram)
   (mass-difference kilogram))
  ((= speed-of-light '(q 299792458 (/ meter second)))
   (= mass-difference (- initial-mass final-mass))
   (= energy (* mass-difference (expt speed-of-light 2)))  ; released
 ))

; ***** light emissions vs. energy band transitions  1998 #FR7
(bohr-model ()
  ()
  ())

; *****
(compton-scatter (frequency)
  ()
  () )

; charged particle moving in a magnetic field
(charge-magnetic (charge field velocity)
  ((field tesla))
  ((= force (* charge (* field velocity)))
   (= radius (/ (* mass velocity) (* charge field)))
 ))

; Thermodynamics

; m = mass, dt = delta temperature, c = specific heat
(heat-transfer (m1 dt1 c1)
  ( (t1f kelvin) (t2f kelvin) (t1 kelvin) (t2 kelvin)
    (dt1 kelvin) (dt2 kelvin) (m1 kilogram) (m2 kilogram)
    (c1 (/ joule (* kilogram degreeK))) (c2 (/ joule (* kilogram degreeK))) )
  ((= dt1 (- t1f t1))
   (= dt2 (- t2f t2))
   (= hm1 (* m1 c1))
   (= hm2 (* m2 c2))
   (= equilibrium-temperature (/ (+ (* hm1 t1) (* hm2 t2)) (+ hm1 hm2)))
   (= energy (* (* m1 c1) dt1))
   (= energy (- (* (* m2 c2) dt2)))
   (= rate (/ energy time)) ))

(ideal-gas (pressure volume)
  ((density (/ kilogram (* meter meter meter)))
   (delta-volume (* meter meter meter))
   (average-speed (/ meter second)) (molar-mass (/ kilogram mole))
   (n-moles mole))
  ((= temperature (/ (* pressure volume)
		     (* n-moles universal-gas-constant)))
   (= universal-gas-constant '(q 8.314472 (/ joule (* mole degree-kelvin))))
   (= n-moles '(q 1 mole))
   (= mass (* n-moles molar-mass))
   (= density (/ mass volume))
   (= work (- (* pressure delta-volume)))  ; assumes pressure const ***
   (= internal-energy (* 3/2 (* n-moles
				(* universal-gas-constant temperature))))
   (= average-speed (sqrt (/ (* 3 (* universal-gas-constant temperature))
			     molar-mass)))
 ))

; heat = heat added to system
; work = work done by system
(thermo-system (initial-energy heat work)
  ((work-absorbed joule) (total-work joule) (energy-change joule)
   (final-energy joule) (initial-energy joule) (heat joule))
  ((= work-absorbed '(q 0 joule))  ; default
   (= total-work (- work work-absorbed))
   (= energy-change (- heat total-work))
   (= final-energy  (+ initial-energy energy-change))
))

(heat-engine (heat work)
  ((efficiency unity) (heat joule) (heat-exhausted joule))
  ((= work (- heat heat-exhausted))
   (= efficiency (/ work heat))
   (= max-efficiency (- 1 (/ temperature-cold temperature-hot)))
  ))

(relativity (c-factor)
  ((relativity-factor unity) (observed-time second)
   (observed-length meter))
  ((= speed-of-light '(q 299792458 (/ meter second)))
   (= velocity (* c-factor speed-of-light))
   (= relativity-factor (sqrt (- 1 (/ (expt velocity 2)
				      (expt speed-of-light 2)))))
   (= observed-time (/ time relativity-factor))
   (= observed-length (* length relativity-factor))
))

; Map coordinates

; 29 Mar 98; 14 Jan 99; 09 Feb 99
(linear-scale (ratio base)
  ()
  ((= width (- xmax xmin))
   (= goalwidth (- goalmax goalmin))
   (= ratio (/ goalwidth width))
   (= base (- goalmin (* ratio xmin)))
   (= y (+ base (* x ratio))) ) )

; 25 Nov 94
(mercator (latitude longitude radius)
  ()
  ((= maplong (* radius longitude))
   (= maplat  (* radius (log (/ (+ 1.0 (sin latitude))
				(cos latitude))))) ) )

(map-scale (xratio xbase yratio ybase)
  ()
  ((= width (- xmax xmin))
   (= height (- ymax ymin))
   (= goalwidth (- goalxmax goalxmin))
   (= goalheight (- goalymax goalymin))
   (= xratio (/ goalwidth width))
   (= yratio (/ goalheight height))
   (= xbase (- goalxmin (* xratio xmin)))
   (= ybase (- goalymin (* yratio ymin)))
   (= goalx (+ xbase (* x xratio)))
   (= goaly (+ ybase (* y yratio))) ) )


; Physical constants, from http://www.physics.nist.gov/
; 2002 CODATA recommended values
; note that the number in parens is +/- error in last digits

; unified atomic mass unit 1.660 538 86(28) x 10-27 kg
; atomic mass unit-electron volt relationship u*c^2  931.494 043(80) x 10+6 eV
; proton mass 1.672 621 71(29) x 10-27 kg
; neutron mass 1.674 927 28(29) x 10-27 kg
; electron mass 9.109 3826(16) x 10-31 kg
; elementary charge  1.602 176 53(14) x 10-19 C
; Avogadro constant 6.022 1415(10) x 10+23 mol-1
; molar gas constant R 8.314 472(15) J mol-1 K-1
; Boltzmann constant k 1.380 6505(24) x 10-23 J K-1
; speed of light in vacuum 299 792 458 m s-1
; Planck constant h 6.626 0693(11) x 10-34 J s
; Planck constant in eV s 4.135 667 43(35) x 10-15 eV s
; hc = 1.9864456023253393E-25 J m     (* 6.6260693E-34 299792458)  * calculated
;    = 1.2398419057638671E+3  eV * nm
;             (/ (* 6.6260693E-34 299792458) 1.60217653E-19)       * calculated
; electric constant E0 = 8.854 187 817... x 10-12 F m-1
;                E0 = permittivity of free space = 1/(4*pi*ke)
; Coulomb's law const ke = 1/(4*pi*e0) = 8.9875517879979115E9 N m2 C-2  *calc
;             (/ 1 (* 4 pi 8.854187817E-12))
; magnetic constant mu0 4pi x 10-7 = 12.566 370 614... x 10-7
; Newtonian constant of gravitation G 6.6742(10) x 10-11 m3 kg-1 s-2
; standard acceleration of gravity 9.806 65 m s-2
; standard atmosphere 101 325 Pa = 101 325 N m-2
; electron volt 1.602 176 53(14) x 10-19 J
; angstrom  1 x 10-10 m
; {220} lattice spacing of silicon 192.015 5965(70) x 10-12 m

; Hardy-Weinberg describes the population distribution of a genetic trait 
; with a single gene locus (ex: Rh factors in human blood types).
; p is the frequency of allele A while q is the frequency of allele a
(population (p)
  ((homozygous-recessive unity) (heterozygous unity)
   (homozygous-dominant unity) (p unity) (q unity))
  ((= 1 (+ p q))
   (= homozygous-recessive (expt q 2))
   (= heterozygous (* 2 (* p q)))
   (= homozygous-dominant (expt p 2)) ))
))

; 01 Apr 04; 05 Apr 04
; find assumed unit of a variable
(defun assumedvartype (var &optional objtype)
  (or (and objtype
	   (cadr (assoc var (vars-units objtype))))
      (assumed-unit var)))

; 29 Jun 06
; problem: assumed units may be expressions
(defun defphyslaw (item)
  (setf (basis-vars (car item)) (second item))
  (setf (vars-units (car item)) (third item))
  (setf (vars-unittypes (car item))
	(mapcar #'(lambda (pair)
		    (list (car pair) (glabstractunit (cadr pair)) ))
		(third item)))
  (setf (equations (car item)) (fourth item)) )

(dolist (item *physlaws*) (defphyslaw item))

(pushnew '(linear-scale linear-scale) *conn-user-laws*)
(pushnew '(map-scale map-scale)  *conn-user-laws*)