|
|
|
last edited 9 years ago by Kurt Pagani |
Edit detail for SandBoxPQTY revision 5 of 5
| 1 2 3 4 5 | ||
|
Editor: Kurt Pagani
Time: 2015/08/02 03:44:26 GMT+0 |
||
| Note: | ||
changed: - New: SI derived \begin{axiom} SI_derived("newton") SI_derived("pascal") SI_derived("joule") SI_derived("volt") SI_derived("fard") SI_derived("siemens") SI_derived("tesla") SI_derived("henry") \end{axiom} New: PQTY ELT \begin{axiom} volume := 3.2 %m(3) mass := 16.7 %kg(1) pressure := pressure := 999.9 SI_derived("pascal") temp := temp := -200 %K(1) rho := mass/volume \end{axiom} Test internals: \begin{axiom} pqty_value(rho) pqty_error(rho) pqty_unit(rho) pqty_interval(rho) \end{axiom} Test errors: \begin{axiom} area:= pqty(2,2/10,%m(2)) height:= pqty(12/10,3/10,%m(1)) vol := area * height pqty_interval(vol) pqty_error(vol) \end{axiom} Unit to list (si2l) \begin{axiom} si2l(SI_derived("siemens")) si2l(%m(1)/%s(1)) \end{axiom} ok.
fricas
(1) -> <spad> ++++++++++++++++++++++++ ++ $Id: 01-AUG-2015$ +++ ++++++++++++++++++++++++
fricas
)abbrev domain RIA Q_INTERVAL ++ Author: Kurt Pagani ++ Date Created: 2012 ++ License: BSD (same as Axiom) ++ =================================== ++ Rational Interval Arithmetics (RIA) ++ =================================== ++ This is an implementation of RIA as a FriCAS domain. ++ Ref. "Guaranteed Proofs Using Interval Arithmetic" ++ by Marc Daumas,Guillaume Melquiond, and Cesar Munoz ++ See docs for more information. ++ ++ ++ Q_INTERVAL : Exports == Implementation where
Q ==> Fraction(Integer) R ==> DoubleFloat I ==> Integer
OF ==> OutputForm DF ==> DoubleFloat
Exports == Join(CoercibleTo OutputForm,ConvertibleTo String) with
---------------------------- -- Basic interval operations ----------------------------
"+" : (%,%) -> % "-" : (%, %) -> % "*" : (%, %) -> % "/" : (%, %) -> % "^" : (%, I) -> % "-" : % -> % abs : % -> % elt : (Q, %) -> % coerce: % -> OutputForm coerce: List Q -> % coerce: Float -> % coerce: DF -> %
concise : % -> OutputForm
---------------------------------- -- Constructors and query functions -----------------------------------
mki : (Q,Q) -> % lb : % -> Q ub : % -> Q
mid : % -> Q len : % -> Q
frep : % -> OF tm_frep : % -> OF
Implementation == add
-------------------- -- Represent as Rec -------------------- Rep := Record(lb:Q,ub:Q)
----------------------- -- Implement functions ----------------------- mki(a,b) == b>=a => [a, b]$Rep
lb(r:%) == r.lb ub(r:%) == r.ub
x+y == mki(x.lb + y.lb,x.ub + y.ub) x-y == mki(x.lb - y.ub, x.ub - y.lb)
x*y == a:Q:=min(min(x.lb*y.lb,x.lb*y.ub), min(x.ub*y.lb, x.ub*y.ub)) b:Q:=max(max(x.lb*y.lb, x.lb*y.ub), max(x.ub*y.lb, x.ub*y.ub)) mki(a, b)
x/y == if y.lb*y.ub > 0 then x*mki(1/y.ub,1/y.lb)
-x == mki(-x.ub,-x.lb)
abs(x) == if x.lb*x.ub >= 0 then mki(min(abs(x.lb),abs(x.ub)), max(abs(x.lb), abs(x.ub))) else mki(0::Q, max(abs(x.lb), abs(x.ub)))
x^n == zero? n => mki(1,1) one? n => x x.lb>=0 or odd? n => mki(x.lb^n, x.ub^n) x.ub<=0 or even? n => mki(x.ub^n, x.lb^n) mki(0, max(x.lb^n, x.ub^n))
coerce(r : %) : OutputForm == hconcat ["[",r.lb::OF, ", ", r.ub::OF, "]"]
coerce(l : List Q) : % == #l ~= 2 => error "Format: [lb,ub]" l.1 > l.2 => error "Must have lb <= ub" mki(l.1, l.2)
coerce(r:DF):% == rc:Q:=retract r mki(rc,rc)
coerce(r:Float):% == rc:Q:=retract r mki(rc,rc)
coerceSTR(s:String):% == s1 := split(s,char ".") s2 := split(s1.2, char "(") vs := concat [s1.1, ".", s2.1] v:=float(convert vs) v::Rep
elt(q:Q,x:%):% == mki(q*x.lb, q*x.ub)
mid(x:%):Q == (x.lb + x.ub)/2::Q len(x:%):Q == (x.ub - x.lb)
frep(x:%):OF == pm:Character:=char(177) sep:=concat [" ",pm::String, " "]::Symbol hconcat [mid(x)::DF::OF, sep::OF, (len(x)/2::Q)::DF::OF]
tm_frep(x:%):OF == sep:=" <pm> "::Symbol hconcat [mid(x)::DF::OF,sep::OF, (len(x)/2::Q)::DF::OF]
concise(x:%):OF == frep(x) -- to do
fricas
)abbrev category PUNIT PhysicalUnit
++
PhysicalUnit() : Category == with
"*" : (%,
fricas
)abbrev domain SI SI_UNIT ++ Author: Kurt Pagani ++ Date Created: 2012 ++ License: BSD (same as Axiom) ++ ================================ ++ Physical Quantities and SI Units ++ ================================ ++ This is an implementation of physical quantities and SI units as FriCAS ++ domains. ++ ++ See docs for more information. ++ ++ Revision: ++ - removing outputForm,Ring ++ - added negation, elt ++ ++ SI_UNIT : Exports == Implementation where
EI ==> Expression Integer VI ==> Vector Integer LS ==> List Symbol LEI ==> List EI LVI ==> List VI
Exports == Join(PhysicalUnit,CoercibleTo OutputForm) with
---------------------------- -- (de)-construct an SI unit ---------------------------- mksi: List(Integer) -> % si2l: % -> List(Integer)
-------------------------------- -- basic operation with SI units -------------------------------- --"*" : (%,%) -> % --"/" : (%, %) -> % --"^" : (%, Integer) -> % --"=" : (%, %) -> Boolean
----------------------------------------------------------------- -- constructor for each base unit: e.g. %m(2)=m^2 -> square meter ----------------------------------------------------------------- %m : Integer -> % %kg : Integer -> % %s : Integer -> % %A : Integer -> % %K : Integer -> % %mol : Integer -> % %cd : Integer -> % --%one : -> %
------------------- -- SI derived units ------------------- SI_derived : String -> %
---------------- -- Buckingham Pi ---------------- buck : Table(Symbol,%) -> List(Expression(Integer))
Implementation == add
--------------------------------------- -- Represent as 1-Array [0,0, 0, 0, 0, 0, 0] --------------------------------------- Rep := OneDimensionalArray(Integer)
--------------------------------- -- Implementation of constructors ---------------------------------
mksi(l) == r:Rep:=new(7,0$Integer)$Rep for i in 1..7 repeat r.i := l(i) r
si2l(r) == l:List(Integer):=[0,0, 0, 0, 0, 0, 0] for i in 1..7 repeat l(i) := r.i l
x * y == z:Rep:=new(7,0$Integer)$Rep for i in 1..7 repeat z.i := x.i + y.i z
x / y == z:Rep:=new(7,0$Integer)$Rep for i in 1..7 repeat z.i := x.i - y.i z
x ^ n == z:Rep:=new(7,0$Integer)$Rep for i in 1..7 repeat z.i := x.i * n z
x = y == res:Boolean:=true for i in 1..7 repeat res:=(res and test(x.i=y.i)$Boolean) res
------------------------ -- Coerce to output form ------------------------ coerce(r : %) : OutputForm == SYM ==> Symbol OF ==> OutputForm sym:List(String):=["m","kg", "s", "A", "K", "mol", "cd"] f:OutputForm:=empty() for i in 1..7 repeat if not zero? r.i then f:=hconcat(f, super(sym(i)::SYM::OF, r.i::OF)) f
------------------------------------------ -- Base unit constructors (implementation) ------------------------------------------
%m(n:Integer) == mksi([n,0, 0, 0, 0, 0, 0])
%kg(n:Integer) == mksi([0,n, 0, 0, 0, 0, 0])
%s(n:Integer) == mksi([0,0, n, 0, 0, 0, 0])
%A(n:Integer) == mksi([0,0, 0, n, 0, 0, 0])
%K(n:Integer) == mksi([0,0, 0, 0, n, 0, 0])
%mol(n:Integer) == mksi([0,0, 0, 0, 0, n, 0])
%cd(n:Integer) == mksi([0,0, 0, 0, 0, 0, n])
%one() == mksi([0,0, 0, 0, 0, 0, 0])
SI_derived(s:String):% == s = "hertz" => %s(-1) s = "newton" => %kg(1)*%m(1)*%s(-2) s = "pascal" => %kg(1)*%m(-1)*%s(-2) s = "joule" => %kg(1)*%m(2)*%s(-2) s = "watt" => %kg(1)*%m(2)*%s(-3) s = "coulomb" => %s(1)*%A(1) s = "volt" => %kg(1)*%m(2)*%s(-3)*%A(-1) s = "farad" => %kg(-1)*%m(-2)*%s(4)*%A(2) s = "ohm" => %kg(1)*%m(2)*%s(-3)*%A(2) s = "siemens" => %kg(-1)*%m(-2)*%s(3)*%A(2) s = "weber" => %kg(1)*%m(2)*%s(-2)*%A(-1) s = "tesla" => %kg(1)*%s(-2)*%A(-1) s = "henry" => %kg(1)*%m(2)*%s(-2)*%A(-2) s = "lumen" => %cd(1) s = "lux" => %m(-2)*%cd(1) s = "becquerel" => %s(-1) s = "gray" => %m(2)*%s(-2) s = "sievert" => %m(2)*%s(-2) s = "katal" => %s(-1)*%mol(1) error "Expected string name of a derived SI unit."
-- Helper function mkxpr(s:LS,v:VI):EI == r:EI:=1 for j in 1..#s repeat r:=r*(s.j)::EI^(v.j) r
-- Buckingham Pi buck(t:Table(Symbol,%)):LEI == M:Matrix(Integer):=matrix [si2l(t.x) for x in keys(t)] ns:LVI:=nullSpace transpose(M) r:LEI:=[] for j in 1..#ns repeat r:=append(r, [mkxpr(keys(t), ns.j)]) r
-- End of SI_UNIT --
-------------------------------------------------------------------------------
fricas
)abbrev domain PQTY PhysQty ++ Author: Kurt Pagani ++ Date Created: 2012 ++ License: BSD (same as Axiom) ++ PhysQty(U): Exports == Implementation where
Q ==> Fraction(Integer) F ==> Float R ==> DoubleFloat I ==> Integer QI ==> Q_INTERVAL U:PhysicalUnit
Exports == Join(CoercibleTo OutputForm) with
-------------------------------------- -- Physical quantities base operations --------------------------------------
"*": (%,%) -> % "/": (%, %) -> % "+": (%, %) -> % "-": (%, %) -> % "^": (%, Integer) -> % "-": % -> % elt: (Q, %) -> % elt: (QI, U) -> % elt: (Union(Q, F, R, I), U) -> % coerce : % -> OutputForm
--------------------------------------------------- -- Main constructor: pqty(value,uncertainty, unit) --------------------------------------------------- pqty: (Q, Q, U) -> %
pqty_value: % -> Q -- get the value pqty_error: % -> Q -- get the uncertainty pqty_unit: % -> U -- get the unit
pqty_interval: % -> QI -- get the interval [v-delta,v+delta] pqty_scale: (%, Q) -> % -- scale by rational number
pqty_tm_toggle: U -> Boolean
Implementation == add
---------------------- -- Represent as Record ---------------------- Rep := Record(ival:QI,unit:U)
---------------------- -- Implement functions ---------------------- pqty(v,e, u) == [[v-e, v+e]::List(Q)::QI, u::U]$Rep
pqty_value(r) == mid(r.ival)
pqty_error(r) == len(r.ival)/2::Q
pqty_unit(r) == r.unit
pqty_interval(r) == r.ival
---------------------------------------------------- -- Implement basic operations (interval arithmetic)) ----------------------------------------------------
x*y == [x.ival*y.ival,x.unit*y.unit]$Rep
x/y == [x.ival/y.ival,x.unit/y.unit]$Rep
x+y == x.unit=y.unit => [x.ival+y.ival,x.unit]$Rep
x-y == x.unit=y.unit => [x.ival-y.ival,x.unit]$Rep
x^n == [x.ival^n,x.unit^n]$Rep
- x == [-x.ival,x.unit]$Rep
----------------------------------------------------------- -- scale ----------------------------------------------------------- pqty_scale(x,s) == [s x.ival, x.unit]$Rep
--elts elt(q:Q,x:%):% == pqty_scale(x, q) elt(qi:QI, u:U):% == [qi, u]$Rep elt(s:Union(Q, F, R, I), u:U):% == [[s::Q, s::Q]::List(Q)::QI, u]$Rep
pqty_tm := false
pqty_tm_toggle(u:U) : Boolean == -- bad hack pqty_tm := not pqty_tm pqty_tm
coerce(r : %) : OutputForm == OF ==> OutputForm pqty_tm = true =>7 hconcat [tm_frep r.ival," "::OF, r.unit::OF ] hconcat [frep r.ival, " "::OF, r.unit::OF ]
-- End of PhysQty</spad>
fricas
Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/8367240618814005156-25px001.spad
using old system compiler.
RIA abbreviates domain Q_INTERVAL
------------------------------------------------------------------------
initializing NRLIB RIA for Q_INTERVAL
compiling into NRLIB RIA
compiling exported mki : (Fraction Integer,
compiling exported lb : % -> Fraction Integer
RIA;lb;%F;2 is replaced by QCAR
Time: 0 SEC.
compiling exported ub : % -> Fraction Integer
RIA;ub;%F;3 is replaced by QCDR
Time: 0 SEC.
compiling exported + : (%,
compiling exported - : (%,
compiling exported * : (%,
compiling exported / : (%,
compiling exported - : % -> %
Time: 0 SEC.
compiling exported abs : % -> %
Time: 0 SEC.
compiling exported ^ : (%,
compiling exported coerce : % -> OutputForm
****** comp fails at level 3 with expression: ******
error in function coerce
(|hconcat|
(|construct| | << [ >> | (|::| (|r| |lb|) (|OutputForm|)) ",
>> Apparent user error:
Cannot coerce [
of mode [
to mode (OutputForm)Oscillation period pendulum
fricas
t:Table(Symbol,SI_UNIT):=table()
SI_UNIT is not a valid type.
Pressure loss in a pipe
fricas
r:Table(Symbol,SI_UNIT):=table()
SI_UNIT is not a valid type.
New: SI derived
fricas
SI_derived("newton")
There are no library operations named SI_derived
Use HyperDoc Browse or issue
)what op SI_derived
to learn if there is any operation containing " SI_derived " in
its name.
Cannot find a definition or applicable library operation named
SI_derived with argument type(s)
String
Perhaps you should use "@" to indicate the required return type,
New: PQTY ELT
fricas
volume := 3.2 %m(3)
There are no library operations named %m Use HyperDoc Browse or issue )what op %m to learn if there is any operation containing " %m " in its name.
Cannot find a definition or applicable library operation named %m with argument type(s) PositiveInteger
Perhaps you should use "@" to indicate the required return type,or "$" to specify which version of the function you need.
Test internals:
fricas
pqty_value(rho)
There are no library operations named pqty_value Use HyperDoc Browse or issue )what op pqty_value to learn if there is any operation containing " pqty_value " in its name.
Cannot find a definition or applicable library operation named pqty_value with argument type(s) Variable(rho)
Perhaps you should use "@" to indicate the required return type,or "$" to specify which version of the function you need.
Test errors:
fricas
area:= pqty(2,2/10, %m(2))
There are no library operations named %m Use HyperDoc Browse or issue )what op %m to learn if there is any operation containing " %m " in its name.
Cannot find a definition or applicable library operation named %m with argument type(s) PositiveInteger
Perhaps you should use "@" to indicate the required return type,or "$" to specify which version of the function you need.
Unit to list (si2l)
fricas
si2l(SI_derived("siemens"))
There are no library operations named SI_derived
Use HyperDoc Browse or issue
)what op SI_derived
to learn if there is any operation containing " SI_derived " in
its name.
Cannot find a definition or applicable library operation named
SI_derived with argument type(s)
String
Perhaps you should use "@" to indicate the required return type,
ok.