Applications, Examples and Libraries

Share your work here

Hi
Two new things recently added to the latest version of Physics available on Maplesoft's R&D Physics webpage are worth mentioning outside the framework of Physics.

  • automaticsimplification. This means that after "Physics:-Setup(automaticsimplification=true)", the output corresponding to every single input (literally) gets automatically simplified in size before being returned to the screen. This is fantastically convenient for interactive work in most situations.

  • Add Physics:-Library:-Assume, to perform the same operations one typically performs with the  assume command, but without the side effect that the variables get redefined. So the variables do not get redefined, they only receive assumptions.

This new Assume implements the concept of an "extended assuming". It permits re-using expressions involving the variables being assumed, expressions that were entered before the assumptions were placed, as well as reusing all the expressions computed while the variables had assumptions, even after removing the variable's assumptions. None of this is possible when placing assumptions using the standard assume. The new routine also permits placing assumptions on global variables that have special meaning, that cannot be redefined, e.g. the cartesian, cylindrical or spherical coordinates sets, or the coordinates of a coordinate spacetime system within the Physics package, etc.

Examples:

 

with(Physics):

This is Physics from today:

Physics:-Version()[2]

`2014, December 9, 16:51 hours`

(1.1)
• 

Automatic simplification is here. At this point automaticsimplification is OFF by default.

Setup(automaticsimplification)

[automaticsimplification = false]

(1.2)

Hence, for instance, if you input the following expression, the computer just echoes your input:

Physics:-`*`(a, c)+Physics:-`*`(a, d)+Physics:-`*`(b, c)+Physics:-`*`(b, d)

a*c+a*d+b*c+b*d

(1.3)

There is however some structure behind (1.3) and, in most situations, it is convenient to have these structures
apparent, in part because they frequently provide hints on how to proceed ahead, but also because a more
compact expression is, roughly speaking, simpler to understand. To see this
automaticsimplification in action,
turn it ON:

Setup(automaticsimplification = true)

[automaticsimplification = true]

(1.4)

Recall this same expression (you could input it with the equation label (1.3) as well) 

Physics:-`*`(a, c)+Physics:-`*`(a, d)+Physics:-`*`(b, c)+Physics:-`*`(b, d)

(c+d)*(a+b)

(1.5)

What happened: this output, as everything else after you set automaticsimplification = true and with no
exceptions, is now further processed with simplify/size before being returned. And enjoy computing with frankly
shorter expressions all around! And no need anymore for "simplify(%, size)" every three or four input lines.

Another  example, typical in computer algebra where expressions become uncomfortably large and difficult to
read: convert the following input to 2D math input mode first, in order to compare what is being entered with the
automatically simplified output on the screen

-Physics:-`*`(Physics:-`*`(Physics:-`*`(3, sin(x)^(1/2)), cos(x)^2), sin(x)^m)+Physics:-`*`(Physics:-`*`(Physics:-`*`(3, sin(x)^(1/2)), cos(x)^2), cos(x)^n)+Physics:-`*`(Physics:-`*`(Physics:-`*`(4, sin(x)^(1/2)), cos(x)^4), sin(x)^m)-Physics:-`*`(Physics:-`*`(Physics:-`*`(4, sin(x)^(1/2)), cos(x)^4), cos(x)^n)

-4*(cos(x)^n-sin(x)^m)*sin(x)^(1/2)*cos(x)^2*(cos(x)^2-3/4)

(1.6)

You can turn automaticsimplification OFF the same way

Setup(automaticsimplification = false)

[automaticsimplification = false]

(1.7)
• 

New Library:-Assume facility; welcome to the world of "extended assuming" :)

 

Consider a generic variable, x. Nothing is known about it

about(x)

x:

  nothing known about this object

 

Each variable has associated a number that depends on the session, and the computer (internally) uses this
number to refer to the variable.

addressof(x)

18446744078082181054

(1.8)

When using the assume  command to place assumptions on a variable, this number, associated to it, changes,
for example:

assume(0 < x and x < Physics:-`*`(Pi, 1/2))

addressof(x)

18446744078179060574

(1.9)

Indeed, the variable x got redefined and renamed, it is not anymore the variable x referenced in (1.8).

about(x)

Originally x, renamed x~:

  is assumed to be: RealRange(Open(0),Open(1/2*Pi))

 


The semantics may seem confusing but that is what happened, you enter x and the computer thinks x~, not x 
anymore.This means two things:

1) all the equations/expressions, entered before placing the assumptions on x using assume, involve a variable x 
that is different than the one that exists after placing the assumptions, and so these previous expressions
cannot
be reused
. They involve a different variable.

2) Also, because, after placing the assumptions using assume, x refers to a different object, programs that depend
on the
x that existed before placing the assumptions will not recognize the new x redefined by assume .

 

For example, if x was part of a coordinate system and the spacetime metric g[mu, nu]depends on it, the new variable x
redefined within assume, being a different symbol, will not be recognized as part of the dependency of "g[mu,nu]." This
posed constant obstacles to working with curved spacetimes that depend on parameters or on coordinates that
have a restricted range. These problems are resolved entirely with this new
Library:-Assume, because it does not
redefine the variables. It only places assumptions on them, and in this sense it works like
assuming , not assume .
As another example, all the
Physics:-Vectors commands look for the cartesian, cylindrical or spherical coordinates
sets
[x, y, z], [rho, phi, z], [r, theta, phi] in order to determine how to proceed, but these variables disappear if you use
assume to place assumptions on them. For that reason, only assuming  was fully compatible with Physics, not assume.

 

To undo assumptions placed using the assume command one reassigns the variable x to itself:

x := 'x'

x

(1.10)

Check the numerical address: it is again equal to (1.8) 

addressof(x)

18446744078082181054

(1.11)

·All these issues get resolved with the new Library:-Assume, that uses all the implementation of the existing 
assume command but with a different approach: the variables being assumed do not get redefined, and hence:
a) you can reuse expressions/equations entered before placing the assumptions, you can also undo the
assumptions and reuse results obtained with assumptions. This is the concept of an
extended assuming. Also,
commands that depend on these assumed variables will all continue to work normally, before, during or after
placing the assumption, because
the variables do not get redefined.

Example:

about(x)

x:

  nothing known about this object

 

So this simplification attempt accomplishes nothing

simplify(arccos(cos(x)))

arccos(cos(x))

(1.12)

Let's assume now that 0 < x and x < (1/2)*Pi

Library:-Assume(0 < x and x < Physics:-`*`(Pi, 1/2))

{x::(RealRange(Open(0), Open((1/2)*Pi)))}

(1.13)

The new command echoes the internal format representing the assumption placed.

a) The address is still the same as (1.8)

addressof(x)

18446744078082181054

(1.14)

So the variable did not get redefined. The system however knows about the assumption - all the machinery of the
assume command is being used

about(x)

Originally x, renamed x:

  is assumed to be: RealRange(Open(0),Open(1/2*Pi))

 


Note that the renaming is to the variable itself - i.e. no renaming.

Hence, expressions entered before placing assumptions can be reused. For example, for (1.12), we now have

simplify(arccos(cos(x)))

x

(1.15)

To clear the assumptions on x, you can use either of Library:-Assume(x=x) or Library:-Assume(clear = {x, ...}) in
the case of many variables being cleared in one go, or in the case of a single variable being cleared:

Library:-Assume(clear = x)

about(x)

x:

  nothing known about this object

 


The implementation includes the additionally functionality, for that purpose add the keyword
additionally 
anywhere in the calling sequence. For example:

Library:-Assume(x::positive)

{x::(RealRange(Open(0), infinity))}

(1.16)

about(x)

Originally x, renamed x:

  is assumed to be: RealRange(Open(0),infinity)

 

Library:-Assume(additionally, x < 1)

{x::(RealRange(Open(0), Open(1)))}

(1.17)

Library:-Assume(x = x)

In summary, the new Library:-Assume command implements the concept of an extended assuming, that can be
turned ON and OFF at will at any moment without changing the variables involved.


Download AutomaticSimplificationAndAssume.mw

 

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft


Partial rectification for the Physics:-Simplify and Physics:-Library:-SortProducts procedures dealing with Fermi annihilation/creation operators

This post will be useful for physicists dealing with Fermi annihilation/creation operators. Physics Package provides plenty of powerful tools for quantum operators handling, however some of them often fail to render correct result.  In particular incorrect behaviour with respect to Fermi annihilation/creation operators is observed for routines Simplify and SortProducts.  In this post I present my procedures S*implifyFermionicOperators and SortProductsFermi which partially solve these issues.

Problems with Physics Package routines

   

Short explanation of custom routines SimplifyFermionicOperators and SortProductsFermi

   

"Details for SimplifyFermionicOperators(z,prefix)"

   

"Details for SortProductsFermi(x,L,prefix)"

   

Weak points

   

Final notes

   


Download FermiCreationAnnihilation.mw

I'd like to pay attention to the article of David Austin "How to Grow and Prune a Classification Tree"

Here is its introduction:

"

It's easy to collect data these days; making sense of it is more work. This article explains a construction in machine learning and data mining called a classification tree. Let's consider an example.

In the late 1970's, researchers at the University of California, San Diego Medical Center performed a study in which they monitored 215 patients following a heart attack. For each patient, 19 variables, such as age and blood pressure, were recorded. Patients were then grouped into two classes, depending on whether or not they survived more than 30 days following the heart attack.

Assuming the patients studied were representative of the more general population of heart attack patients, the researchers aimed to distill all this data into a simple test to identify new patients at risk of dying within 30 days of a heart attack.

By applying the algorithm described here, Breiman, Freidman, Olshen, and Stone, created a test consisting of only three questions---what is the patient's minimum systolic blood pressure within 24 hours of being admitted to the hospital, what is the patient's age, does the patient exhibit sinus tachycardia---to identify patients at risk. In spite of its simplicity, this test proved to be more accurate than any other known test. In addition, the importance of these three questions indicate that, among all 19 variables, these three factors play an important role in determining a patient's chance of surviving.

Besides medicine, these ideas are applicable to a wide range of problems, such as identifying which loan applicants are likely to default and which voters are likely to vote for a particular political party.

In what follows, we will describe the work of Breiman and his colleagues as set out in their seminal book Classification and Regression Trees. Theirs is a very rich story, and we will concentrate on only the essential ideas"

It would be interesting to compare this approach with discriminant analysis. Hope somebody of  the Maple developers will give a concrete example on this theme with Maple.

The procedure  NumbersGame  generalizes the well-known 24 game  (implementation in Maple see here), as well as related issues (see here and here).

 

Required parameters of the procedure:

Result is an integer or a fraction of any sign.

Numbers is a list of positive integers.

 

Optional parameters:

Operators is a list of permitted arithmetic operations. By default  Operators is  ["+","-","*","/"]

NumbersOrder is a string. It is equal to "strict order" or "arbitrary order" . By default  NumbersOrder is "strict order"

Parentheses is a symbol  no  or  yes . By default  Parentheses is  no 

 

The procedure puts the signs of operations from the list  Operators  between the numbers from  Numbers  so that the result is equal to Result. The procedure finds all possible solutions. The global  M  saves the list of the all solutions.

 

Code the procedure:

restart;

NumbersGame:=proc(Result::{integer,fraction}, Numbers::list(posint), Operators::list:=["+","-","*","/"], NumbersOrder::string:="strict order", Parentheses::symbol:=no)

local MyHandler, It, K, i, P, S, n, s, L, c;

global M;

uses StringTools, ListTools, combinat; 

 MyHandler := proc(operator,operands,default_value)

      NumericStatus( division_by_zero = false );

      return infinity;

   end proc;

   NumericEventHandler(division_by_zero=MyHandler); 

if Parentheses=yes then 

It:=proc(L1,L2)

local i, j, L;

for i in L1 do

for j in L2 do

L[i,j]:=seq(Substitute(Substitute(Substitute("( i Op j )","i",convert(i,string)),"j",convert(j,string)),"Op",Operators[k]), k=1..nops(Operators));

od; od;

L:=convert(L, list);

end proc; 

P:=proc(L::list)

local n, K, i, M1, M2, S;

n:=nops(L);

if n=1 then return L else

for i to n-1 do

M1:=P(L[1..i]); M2:=P(L[i+1..n]);

K[i]:=seq(seq(It(M1[j], M2[k]), k=1..nops(M2)), j=1..nops(M1))

od; fi;

K:=convert(K,list);

end proc;

if NumbersOrder="arbitrary order" then S:=permute(Numbers); K:=[seq(op(Flatten([op(P(s))])), s=S)] else  K:=[op(Flatten([op(P(Numbers))]))] fi; 

else 

if NumbersOrder="strict order" then

K:=[convert(Numbers[1],string)];

for i in Numbers[2..-1] do

K:=[seq(seq(cat(k, Substitute(Substitute(" j i","j",convert(j,string)),"i",convert(i,string))), k in K), j in Operators)]

od;   else 

S:=permute(Numbers);

for s in S do

L:=[convert(s[1],string)];

for i in s[2..-1] do

L:=[seq(seq(cat(k, Substitute(Substitute(" j i","j",convert(j,string)),"i",convert(i,string))), k in L), j in Operators)]

od; K[s]:=op(L) od; K:=convert(K,list) fi;  

fi; 

M:='M'; c:=0;

for i in K do

if parse(i)=Result then c:=c+1; if Parentheses=yes then M[i]:= convert(SubString(i,2..length(i)-1),symbol)=convert(Result,symbol) else M[i]:=convert(i,symbol)=convert(Result,symbol) fi; fi;

od; 

if c=0 then M:=[]; return `No solutions` else M:=convert(M,list);  op(M) fi; 

end proc:

 

Examples of use.

 

Example 1:

NumbersGame(1/20, [$ 1..9]);

     1 * 2 - 3 + 4 / 5 / 6 * 7 / 8 * 9 = 1/20

 

Example 2. Numbers in the list  Numbers  may be repeated and permitted operations can be reduced:

NumbersGame(15, [3,3,5,5,5], ["+","-"]);

         3 - 3 + 5 + 5 + 5 = 15

 

Example 3. 

NumbersGame(10, [1,2,3,4,5]);

         1 + 2 + 3 * 4 - 5 = 10

If the order of the number in Numbers is arbitrary, then the number of solutions is greatly increased (10 solutions displayed):

NumbersGame(10, [1,2,3,4,5], "arbitrary order"):

nops(M);

for i to 10 do

M[1+50*(i-1)] od;

If you use the parentheses, the number of solutions will increase significantly more (10 solutions displayed):

NumbersGame(10, [1,2,3,4,5], "arbitrary order", yes):

nops(M);

for i to 10 do

M[1+600*(i-1)] od;

Game.mws

 

 

 

We received an interesting and timely submission to the Maple Application Center this morning that I think people might be interested in.  It's called:

The Comet 67P/Churyumov-Gerasimenko, Rosetta & Philae, by Dr. Ahmed Baroudy. From the abstract:

Our plan is rather a modest one since all we want is to get , by calculations, specific data concerning the comet and its lander.
We shall take a simplified model and consider the comet as a perfect solid sphere to which we can apply Newton's laws.

We want to find:

I- the acceleration on the comet surface ,
II- its radius,
III- its density,
IV- the velocity of Philae just after the 1st bounce off the comet (it has bounced twice),
V- the time for Philae to reach altitude of 1000 m above the comet .

We shall compare our findings with the already known data to see how close our simplified mathematical model findings are to the duck-shaped comet already known results.
It turned out that our calculations for a sphere shaped comet are very close to the already known data.

Click on the link above if you want to take a look.

 

eithne

What is Groebner? That was asked in different forms several times in MaplePrimes and MathStackExchange (for example, see http://math.stackexchange.com/questions/3550/using-gr?bner-bases-for-solving-polynomial-equations ). In view of this I think the presented post on Groebner basis will be useful. This post consists of two parts: its mathematical background and examples of solutions of polynomial systems by hand and with Maple.

Let us start. Up to Wiki http://en.wikipedia.org/wiki/Gr%C3%B6bner_basis ,Groebner basis computation can be seen as a multivariate, non-linear generalization of both Euclid's algorithm for computing polynomial greatest common divisors, and Gaussian elimination for linear systems. This is implemented in Maple trough the Groebner package.
The simplest introduction to the topic I know is a well-written book of Ivan Arzhantsev (https://zbmath.org/?q=an:05864974) which includes the proofs of all the claimed theorems and the solutions of all the exercises. Here is its digest groebner.pdf done by me (The reader is assumed to be familiar with the ideal notion and ring notion (one may refresh her/his knowledge, looking in http://en.wikipedia.org/wiki/Ideal_%28ring_theory%29)). It should be noted that there is no easy reading about this serious matter.
Referring to the digest as appropriate, we solve the system
S:={a*b = c^2+c, a^2 =a+ b*c, a*c = b^2+b} by hand and with the Groebner package.
For the order a > b > c we construct its ideal
J(S):=<f1 = a*b-c^2-c,f2 = a^2-a-b*c, f3 = a*c-b^2-b>.
The link between f1 and f2 gives
f1*a-f2*b = (-c^2-c)*a + (a + b*c)*b = a*b -a*c + b^2*c -
a*c^2 =f4.
The reduction with f1 produces
f4 ->-a*c^2- a*c + b^2*c + c^2 +c =: f4.
Now the reduction with f3 produces
f4 -> -b^2- b - b*c +c^2 + c =:f4.
The link between f2 and f3 gives:
f2*c - f3*a = a*b +a*b^2 -a*c -b*c^2 = f5.
The reduction with f1 produces
f5 -> -a*c + c*b +c^2 +c =:f5.
The reduction with f3 produces
f5 -> -b^2 -b + c*b +c^2 +c =:f5.
The reduction with f4 produces
f5 -> 2b*c =: f5.
The link between f1 and f3
f1*c - f3*b = b^3 + b^2 -c^3 -c^2=:f6.
The reduction with f4 produces
f6 -> 2b*c + 2b*c^2 -2c^3 -2c^2=:f6.
At last, we reduce f6 by f5, obtaining f6:= -2c^3 -2c^2.
We see the minimal reduced Groebner basis of S consists of
a^2 -a -b*c, -b^2 -b- b*c +c^2 +c, -2c^3 - 2c^2.
Now we find the solution set of the system under consideration. The equation -2c^3 - 2c^2 = 0 implies
c=0, c=0, c=-1. The the equation -b^2 - b - b*c +c^2 + c = 0 gives
b = 0 , b = -1, b = 0, b = -1, b = 0, b = 0 respectively.
At last, knowing b and c, we find a from a^2 -a -b*c = 0.
Hence,
[{a = 0, b = 0, c = 0}, {a = 1, b = 0, c = 0}, {a = 0, b = -1, c = 0}], [{a = 0, b = 0, c = 0}, {a = 1, b = 0, c = 0}, {a = 0, b = -1, c = 0}], [{a = 0, b = 0, c = -1}].
The solution of the system under consideration by the Groebner package is somewhat different because Maple does not find the minimal reduced Groebner basis directly.

 

with(Groebner):

[b*c, a*c-c^2-c, b^2-c^2+b-c, a*b-c^2-c, a^2-a, c^3+c^2]

(1)

GB2 := remove(has, GB1, a);

[b*c, b^2-c^2+b-c, c^3+c^2]

(2)

GB3 := Basis(GB2, lexdeg([b, c]))

[b*c, b^2-c^2+b-c, c^3+c^2]

(3)

op(remove(has, GB3, {b}))

c^3+c^2

(4)

solc := solve(c^3+c^2);

-1, 0, 0

(5)

solb := [seq(op(map(`union`, [solve(eval(GB3, c = i), {b})], {c = i})), i = solc)]

[{b = 0, c = -1}, {b = -1, c = 0}, {b = 0, c = 0}, {b = -1, c = 0}, {b = 0, c = 0}]

(6)

sol := seq(op(map(`union`, [solve(eval(GB1, i))], i)), i = solb)

{a = 0, b = 0, c = -1}, {a = 0, b = -1, c = 0}, {a = 0, b = 0, c = 0}, {a = 1, b = 0, c = 0}, {a = 0, b = -1, c = 0}, {a = 0, b = 0, c = 0}, {a = 1, b = 0, c = 0}

(7)

NULL

S2 := {a*c-b^2-b, a*b-c^2-c, a^2-b*c+a}:

GB1 := Basis(S2, lexdeg([a, b, c]))

[a*c+b*c+c^2+c, b^2+b*c+c^2+b+c, a*b-c^2-c, a^2-b*c+a]

(8)

GB2 := remove(has, GB1, a)

[b^2+b*c+c^2+b+c]

(9)

sol1 := solve(GB2, b)

{b = -(1/2)*c-1/2+(1/2)*(-3*c^2-2*c+1)^(1/2)}, {b = -(1/2)*c-1/2-(1/2)*(-3*c^2-2*c+1)^(1/2)}

(10)

map(proc (c) options operator, arrow; `union`(c, sol1[1]) end proc, map(proc (C) options operator, arrow; solve(C, {a}) end proc, eval(S2, sol1[1])))

{{a = 2*c*(c+1)/(-c-1+(-3*c^2-2*c+1)^(1/2)), b = -(1/2)*c-1/2+(1/2)*(-3*c^2-2*c+1)^(1/2)}, {a = -1/2-(1/2)*(1+2*c*(-3*c^2-2*c+1)^(1/2)-2*c^2-2*c)^(1/2), b = -(1/2)*c-1/2+(1/2)*(-3*c^2-2*c+1)^(1/2)}, {a = -1/2+(1/2)*(1+2*c*(-3*c^2-2*c+1)^(1/2)-2*c^2-2*c)^(1/2), b = -(1/2)*c-1/2+(1/2)*(-3*c^2-2*c+1)^(1/2)}, {a = -(1/2)*c-1/2-(1/2)*(-3*c^2-2*c+1)^(1/2), b = -(1/2)*c-1/2+(1/2)*(-3*c^2-2*c+1)^(1/2)}}

(11)

``

 

Download groebner.mw

I'd like to pay attention to the recent article "The Misfortunes of a Trio of Mathematicians Using Computer Algebra Systems. Can We Trust in Them?"

In particular, the authors consider the integral

int(abs(exp(2*Pi*Ix)+exp(2*Pi*I*y)),[x=0..1,y=0..1]),

stating "Both Mathematica and Maple return zero as the answer to this calculation. Yet this cannot be correct, because the integrand is clearly positive and nonzero in the indicated region". Unfortunately, they give only the Mathematica command to this end.

Of course, the integral under consideration is complicated so the the simple-minded trials

int(evalc(abs(exp((2*Pi*I)*x)+exp((2*Pi*I)*y))), [x = 0 .. 1, y = 0 .. 1]);

and

VectorCalculus:-int(evalc(abs(exp((2*Pi*I)*x)+exp((2*Pi*I)*y))), [x,y]=Rectangle( 0 .. 1, 0 .. 1));

fail. However,this can be found with Maple (I think with Mathematica too.) in such a way.

 

A := evalc(abs(exp((2*Pi*I)*x)+exp((2*Pi*I)*y)))

((cos(2*Pi*x)+cos(2*Pi*y))^2+(sin(2*Pi*x)+sin(2*Pi*y))^2)^(1/2)

(1)

NULL

B := simplify(A, trig)

(2*cos(2*Pi*x)*cos(2*Pi*y)+2+2*sin(2*Pi*x)*sin(2*Pi*y))^(1/2)

(2)

op(B)[1]

2*cos(2*Pi*x)*cos(2*Pi*y)+2+2*sin(2*Pi*x)*sin(2*Pi*y)

(3)

combine(op(B)[1], x)

2*cos(2*Pi*x-2*Pi*y)+2

(4)

C := eval(B, op(B)[1] = combine(op(B)[1], x))

(2*cos(2*Pi*x-2*Pi*y)+2)^(1/2)

(5)

int(C, [x = 0 .. 1, y = 0 .. 1])

4/Pi

(6)

``

 

Download int.mw

 

 

Hi, we recently put together a web video on how memes spread on the internet using several visualizations generated from Maple 18:

http://youtu.be/vEhAkEPwESI

Found the new ability to specify a background image for plots to be very helpful.

We have just released a new version of MapleSim.

MapleSim 7 makes it substantially easier to explore and validate designs, create and manage libraries of custom components, and use your MapleSim models with other tools. It includes:

  • Easy model investigation. A new Results Manager gives you greater flexibility when it comes to investigating your simulation results, including the ability to compare simulation runs on the same axes, instantly plot both probed and unprobed variables, and easily create custom plots.
  • Convenient library creation. With MapleSim 7, it is significantly easier to create, manage, and share libraries of custom components.
  • Improved Modelica support. MapleSim 7 expands the support of the Modelica language so that more Modelica definitions can be used directly inside MapleSim.

We have also updated and expanded the MapleSim 7 family of add-on products:

  • The new MapleSim Battery Library, which is available as a separate add-on, allows you to incorporate physics-based predictive models of battery cells into your system models so you can take battery behavior into account early in the design process. 
  • The MapleSim Connector for FMI, which allows engineers to share very efficient, high-fidelity models created in MapleSim with other modeling tools, has been expanded to support more export formats for co-simulation and model exchange.

See What’s New in MapleSim 7 for more information about these and other improvements in MapleSim.

 

eithne

Last week the Physics package was presented in a talk at the Perimeter Institute for Theoretical Physics and in a combined Applied Mathematics and Physics Seminar at the University of Waterloo. The presentation at the Perimeter Institute got recorded. It was a nice opportunity to surprise people with the recent advances in the package. It follows the presentation with sections closed, and at the end there is a link to a pdf with the sections open and to the related worksheet, used to run the computations in real time during the presentation.

COMPUTER ALGEBRA FOR THEORETICAL PHYSICS

 

  

Generally speaking, physicists still experience that computing with paper and pencil is in most cases simpler than computing on a Computer Algebra worksheet. On the other hand, recent developments in the Maple system implemented most of the mathematical objects and mathematics used in theoretical physics computations, and dramatically approximated the notation used in the computer to the one used in paper and pencil, diminishing the learning gap and computer-syntax distraction to a strict minimum. In connection, in this talk the Physics project at Maplesoft is presented and the resulting Physics package illustrated tackling problems in classical and quantum mechanics, general relativity and field theory. In addition to the 10 a.m lecture, there will be a hands-on workshop at 1pm in the Alice Room.

 

... Why computers?

 

 

We can concentrate more on the ideas instead of on the algebraic manipulations

 

We can extend results with ease

 

We can explore the mathematics surrounding a problem

 

We can share results in a reproducible way

 

Representation issues that were preventing the use of computer algebra in Physics

 

 

Notation and related mathematical methods that were missing:


coordinate free representations for vectors and vectorial differential operators,

covariant tensors distinguished from contravariant tensors,

functional differentiation, relativity differential operators and sum rule for tensor contracted (repeated) indices

Bras, Kets, projectors and all related to Dirac's notation in Quantum Mechanics

 

Inert representations of operations, mathematical functions, and related typesetting were missing:

 

inert versus active representations for mathematical operations

ability to move from inert to active representations of computations and viceversa as necessary

hand-like style for entering computations and texbook-like notation for displaying results

 

Key elements of the computational domain of theoretical physics were missing:

 

ability to handle products and derivatives involving commutative, anticommutative and noncommutative variables and functions

ability to perform computations taking into account custom-defined algebra rules of different kinds

(problem related commutator, anticommutator, bracket, etc. rules)

Vector and tensor notation in mechanics, electrodynamics and relativity

   

Dirac's notation in quantum mechanics

   

 

• 

Computer algebra systems were not originally designed to work with this compact notation, having attached so dense mathematical contents, active and inert representations of operations, not commutative and customizable algebraic computational domain, and the related mathematical methods, all this typically present in computations in theoretical physics.

• 

This situation has changed. The notation and related mathematical methods are now implemented.

 

Tackling examples with the Physics package

 

Classical Mechanics

 

Inertia tensor for a triatomic molecule

 

 

Problem: Determine the Inertia tensor of a triatomic molecule that has the form of an isosceles triangle with two masses m[1] in the extremes of the base and mass m[2] in the third vertex. The distance between the two masses m[1] is equal to a, and the height of the triangle is equal to h.

Solution

   

Quantum mechanics

 

Quantization of the energy of a particle in a magnetic field

 


Show that the energy of a particle in a constant magnetic field oriented along the z axis can be written as

H = `&hbar;`*`&omega;__c`*(`#msup(mi("a",mathcolor = "olive"),mo("&dagger;"))`*a+1/2)

where `#msup(mi("a",mathcolor = "olive"),mo("&dagger;"))`and a are creation and anihilation operators.

Solution

   

The quantum operator components of `#mover(mi("L",mathcolor = "olive"),mo("&rarr;",fontstyle = "italic"))` satisfy "[L[j],L[k]][-]=i `&epsilon;`[j,k,m] L[m]"

   

Unitary Operators in Quantum Mechanics

 

(with Pascal Szriftgiser, from Laboratoire PhLAM, Université Lille 1, France)

A linear operator U is unitary if 1/U = `#msup(mi("U"),mo("&dagger;"))`, in which case, U*`#msup(mi("U"),mo("&dagger;"))` = U*`#msup(mi("U"),mo("&dagger;"))` and U*`#msup(mi("U"),mo("&dagger;"))` = 1.Unitary operators are used to change the basis inside an Hilbert space, which physically means changing the point of view of the considered problem, but not the underlying physics. Examples: translations, rotations and the parity operator.

1) Eigenvalues of an unitary operator and exponential of Hermitian operators

   

2) Properties of unitary operators

   

3) Schrödinger equation and unitary transform

   

4) Translation operators

   

Classical Field Theory

 

The field equations for a quantum system of identical particles

 

 

Problem: derive the field equation describing the ground state of a quantum system of identical particles (bosons), that is, the Gross-Pitaevskii equation (GPE). This equation is particularly useful to describe Bose-Einstein condensates (BEC).

Solution

   

The field equations for the lambda*Phi^4 model

   

Maxwell equations departing from the 4-dimensional Action for Electrodynamics

   

General Relativity

 

Given the spacetime metric,

g[mu, nu] = (Matrix(4, 4, {(1, 1) = -exp(lambda(r)), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -r^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -r^2*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = exp(nu(r))}))

a) Compute the trace of

"Z[alpha]^(beta)=Phi R[alpha]^(beta)+`&Dscr;`[alpha]`&Dscr;`[]^(beta) Phi+T[alpha]^(beta)"

where `&equiv;`(Phi, Phi(r)) is some function of the radial coordinate, R[alpha, `~beta`] is the Ricci tensor, `&Dscr;`[alpha] is the covariant derivative operator and T[alpha, `~beta`] is the stress-energy tensor

T[alpha, beta] = (Matrix(4, 4, {(1, 1) = 8*exp(lambda(r))*Pi, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 8*r^2*Pi, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 8*r^2*sin(theta)^2*Pi, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 8*exp(nu(r))*Pi*epsilon}))

b) Compute the components of "W[alpha]^(beta)"" &equiv;"the traceless part of  "Z[alpha]^(beta)" of item a)

c) Compute an exact solution to the nonlinear system of differential equations conformed by the components of  "W[alpha]^(beta)" obtained in b)

Background: paper from February/2013, "Withholding Potentials, Absence of Ghosts and Relationship between Minimal Dilatonic Gravity and f(R) Theories", by P. Fiziev.

a) The trace of "  Z[alpha]^(beta)=Phi R[alpha]^(beta)+`&Dscr;`[alpha]`&Dscr;`[]^(beta) Phi+T[alpha]^(beta)"

   

b) The components of "W[alpha]^(beta)"" &equiv;"the traceless part of " Z[alpha]^(beta)"

   

c) An exact solution for the nonlinear system of differential equations conformed by the components of  "W[alpha]^(beta)"

   

The Physics Project

 

 

"Physics" is a software project at Maplesoft that started in 2006. The idea is to develop a computational symbolic/numeric environment specifically for Physics, targeting educational and research needs in equal footing, and resembling as much as possible the flexible style of computations used with paper and pencil. The main reference for the project is the Landau and Lifshitz Course of Theoretical Physics.

 

A first version of "Physics" with basic functionality appeared in 2007. Since then the package has been growing every year, including now, among other things, a searcheable database of solutions to Einstein equations and a new dedicated programming language for Physics.

 

Since August/2013, weekly updates of the Physics package are distributed on the web, including the new developments related to our plan as well as related to people's feedback.

 

 

Presentation_at_PI_and_UW.pdf     Presentation_at_PI_and_UW.mw

 

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

The Embedded Components are containers that currently use industries for modeling complex systems to find viable solutions in real time and thus avoid huge wait times and overload our computer; by this paper should show you how to implement a dynamic worksheet through Embedded Components in Maple; it goes from finding solutions to ordinary differential equations partial; which interact with the researcher using different parameters.
Using graphical programming will find immediate solutions to selected problems in science and engineering criteria of variability and boundary conditions evolving development with buttons on multiple actions.

 

cimac_2014.pdf

(in spanish)

Solutions_of_Differential_Equations_with_Embedded_Components.mw

 

Lenin Araujo Castillo

Physics Pure

Computer Science

 

Maplesoft regularly hosts live webinars on a variety of topics. Below you will find details on some upcoming webinars we think may be of interest to the MaplePrimes community.  For the complete list of upcoming webinars, visit our website.

Maplesoft Solutions for Math Education

This webinar will demonstrate how Maplesoft’s solutions for mathematics education help teachers bring complex problems to life, allow students to focus on concepts rather than the mechanics of solutions, and offer students the necessary practice to master the concepts being taught.

Key takeaways include:

• How to quickly and painlessly place incoming students in the correct math courses

• How you can use hundreds of intuitive Clickable Math tools to demonstrate and explore up to advanced-level problems and algorithms in the classroom

• How to automate your testing and assessment needs, specifically for math courses

• How to bring your STEM courses to life in an online environment

To join us for the live presentation, please click here to register.

Introduction to Maple T.A. Placement Test Suite 10

This webinar will provide an overview and demonstration of the latest release of the Maple T.A. MAA Placement Test Suite. A result of the ongoing partnership between the Mathematical Association of America (MAA) and Maplesoft, this product gives you the ability to provide the renowned MAA placement tests in an online testing environment. Learn how the Maple T.A. MAA Placement Test Suite can greatly simplify your placement process and explore the latest additions, including a streamlined interface and new tests to determine your students’ readiness for Precalculus and Algebra courses.

To join us for the live presentation, please click here to register.

There is also a recording available from another live webinar we did earlier this month: Introduction to Maple T.A. 10.

This application calculates the number of photons reaching a camera sensor for a given exposure. A blackbody model of the sun is generated. The "Sunny 16" rule for exposure is demonstrated. Calculations are done using units.Photon_Exposure_Array.mw

Photon ExposureNULLNULL

Blackbody Model of the Sun

    h := Units:-Standard:-`*`(Units:-Standard:-`*`(0.6626069e-33, Units:-Standard:-`^`(Unit('m'), 2)), Units:-Standard:-`*`(Unit('kg'), Units:-Standard:-`/`(Unit('s')))): 

Plank Constant       

  kb := Units:-Standard:-`*`(Units:-Standard:-`*`(0.1380650e-22, Units:-Standard:-`*`(Units:-Standard:-`^`(Unit('m'), 2), Units:-Standard:-`/`(Units:-Standard:-`^`(Unit('s'), 2)))), Units:-Standard:-`*`(Unit('kg'), Units:-Standard:-`/`(Unit('K')))): 

Boltzman Constant  

c := Units:-Standard:-`*`(0.2997925e9, Units:-Standard:-`*`(Unit('m'), Units:-Standard:-`/`(Unit('s')))):  ``

Light Speed

Rsun := Units:-Standard:-`*`(Units:-Standard:-`*`(6.955, Units:-Standard:-`^`(10, 8)), Unit('m')): ``

Sun Radius  

Re_orb := Units:-Standard:-`*`(Units:-Standard:-`*`(1.496, Units:-Standard:-`^`(10, 11)), Unit('m')): ``

Earth Orbit

Tsun := Units:-Standard:-`*`(5800, Unit('K')): ``

Sun Color Temperature     

 tf_atm := .718: 

Transmission Factor  

 

Sun: Spectral Radiant Exitance to Earth: Spectral Irradiance                   

  "M(lambda):=(2*Pi*h*c^(2))/((lambda)^(5))*1/((e)^((h*c)/(lambda*kb*Tsun))-1)*(Rsun/(Re_orb))^(2)*tf_atm:" NULL

evalf(M(Units:-Standard:-`*`(555, Unit('nm')))) = 1277414308.*Units:-Unit(('kg')/(('m')*('s')^3))"(->)"1.277414308*Units:-Unit(('W')/(('nm')*('m')^2))NULL

Photopic Relative Response VP vs λ

 

csvFile := FileTools[Filename]("/VPhotopic.csv")NULL = "VPhotopic.csv"NULL

VPdata := ImportMatrix(csvFile) = Vector(4, {(1) = ` 471 x 2 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})NULLNULL

 

`&lambda;P` := [seq(1 .. 4000)]:

VP := ArrayInterpolation(VPdata, `&lambda;P`):             (ArrayInterpolation for x,y data VPdata returns y' for new x data lambdaP)

NULLVParray := [`$`([`&lambda;P`[n], VP[n]], n = 1 .. 4000)]:                     

Mearth := [`$`([n, Units:-Standard:-`*`(Units:-Standard:-`*`(M(Units:-Standard:-`*`(n, Unit('nm'))), Unit('nm')), Units:-Standard:-`*`(Units:-Standard:-`^`(Unit('s'), 3), Units:-Standard:-`/`(Unit('kg'))))], n = 1 .. 4000)]:````

``

dualaxisplot(plot([Mearth], lambda = 300 .. 900, style = line, color = [blue], labels = ["&lambda; (nm)", "M (W/nm m^2)"], title = "Spectral Radiant Exitance of the Sun", titlefont = ["ARIAL", 15], legend = [Exitance], size = [800, 300]), plot([VParray], style = line, color = [green], labels = ["&lambda; (nm)", "Relative Response"], legend = [Units:-Standard:-`*`(Units:-Standard:-`*`(Photopic, Relative), Response)]))

 

``

 

 

 

Illuminance in Radiometric and Photometric Units:

E__r := sum(Units:-Standard:-`*`(M(Units:-Standard:-`*`(lambda, Unit('nm'))), Unit('nm')), lambda = 200 .. 4000) = 984.7275549*Units:-Unit(('kg')/('s')^3)"(->)"984.7275549*Units:-Unit(('W')/('m')^2)NULL

NULL

E__po := Units:-Standard:-`*`(Units:-Standard:-`*`(683.002, Units:-Standard:-`*`(Unit('lm'), Units:-Standard:-`/`(Unit('W')))), sum(Units:-Standard:-`*`(Units:-Standard:-`*`(VP[lambda], M(Units:-Standard:-`*`(lambda, Unit('nm')))), Unit('nm')), lambda = 200 .. 4000)) = HFloat(91873.47376063903)*Units:-Unit('lx')NULL

Translation from Illuminance to Luminance for Reflected Light;

 

Object Reflectance          R__o:      

Object Luminance           L__po := proc (R__o) options operator, arrow; R__o*E__po/(Pi*Unit('sr')) end proc:                evalf(L__po(1)) = HFloat(29244.234968360346)*Units:-Unit(('cd')/('m')^2) 

 

Illuminance of a Camera Sensor  Eps applied for time texp determines Luminous Exposure Hp;

Ideal Illuminance is determined by the exposure time texp, effective f-number N and to a less extent the angle to the optical axis θ;

 

• 

H       Luminous Exposure

• 

Eps     Illuminance to the Camera

• 

N                                               Effective F-Number

• 

texp             Exposure Time

• 

θ        Angle to the Optical Axis    

 

E__ps_ideal = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), L__po), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2)))):

H__p_ideal = Units:-Standard:-`*`(E__ps_ideal, t__exp):

 

The camera meter determines the exposure time texp to balance the object luminance, reflectance and effective f-number. It does this based on an internal constant k and the camera ISO s.

• 

s        ISO Gain (Based on saturation at 3 stops above the average scene luminance)

• 

k       Reflected Light Meter Calibration Constant      k__m := Units:-Standard:-`*`(Units:-Standard:-`*`(12.5, Unit('lx')), Unit('s')):  

                                                                                                  for Nikon, Canon and Sekonic

• 

c        Incident Light Meter Calibration Constant       c__m := Units:-Standard:-`*`(Units:-Standard:-`*`(250, Unit('lx')), Unit('s')):        

                                                                                                  for Sekonic with flat domeNULL

N^2/t__exp = `#mrow(mi("\`E__po\`"),mo("&sdot;"),mi("s"))`/c__m                        (Incident Light Meter)  NULL 

Units:-Standard:-`*`(Units:-Standard:-`^`(N, 2), Units:-Standard:-`/`(t__exp)) = Units:-Standard:-`*`(`#mrow(mi("\`L__po\`"),mo("&sdot;"),mi("s"))`, Units:-Standard:-`/`(k__m)):                        (Reflected Light Meter)

NULL

Solve for H in terms of the Camera Meter Constant k and s

 

Es = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Lo), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2)))): NULL

t = Units:-Standard:-`*`(Units:-Standard:-`*`(km, Units:-Standard:-`^`(N, 2)), Units:-Standard:-`/`(Units:-Standard:-`*`(Lo, s))):NULL

NULL

NULL

H = Es*t

H = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Lo), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2)))), Units:-Standard:-`*`(Units:-Standard:-`*`(km, Units:-Standard:-`^`(N, 2)), Units:-Standard:-`/`(Units:-Standard:-`*`(Lo, s))))"(=)"H = (1/4)*Pi*cos(theta)^4*km/sNULLNULL

 t = H/Es

t = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), km), Units:-Standard:-`/`(s))), Units:-Standard:-`/`(Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Lo), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2))))))"(=)"t = km*N^2/(Lo*s)NULLNULL

H__p := proc (s, theta) options operator, arrow; (1/4)*Pi*k__m*cos(theta)^4/s end proc:                                              

  evalf(H__p(100, 0)) = 0.9817477044e-1*Units:-Unit(('cd')*('s')/('m')('radius')^2)"(->)"0.9817477044e-1*Units:-Unit(('lx')*('s'))NULL

 

Note:  Meters are typically set for a scene reflectance 3 stops below 100% or 12.5%.

           

  E__ps := proc (N, R__o, theta) options operator, arrow; (1/4)*Pi*Unit('sr')*R__o*E__po*cos(theta)^4/(Pi*Unit('sr')*N^2) end proc:               

 evalf(E__ps(16, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)), 0)) = HFloat(11.215023652421756)*Units:-Unit('lx')                                                                                                   

t__exp_ideal := proc (N, s, R__o) options operator, arrow; H__p(s, theta)/E__ps(N, R__o, theta) end proc:                                     

  evalf(t__exp_ideal(16, 100, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)))) = HFloat(0.008753862094289947)*Units:-Unit('s') NULL NULL

 

 

Actual exposure time includes typical lens losses;

 m := Units:-Standard:-`/`(80):``

Magnification  

  T := .9:``

Lens Transmittance

 F := 1.03:``

Lens Flare

V := 1: ``

Vignetting

 

                                                  ``

Total Lens Efficiency

q := Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(T, F), V), Units:-Standard:-`^`(Units:-Standard:-`+`(1, Units:-Standard:-`-`(m)), 2)):                                      evalf(q) = .9039698438NULL

 

Replacing Eps with q*Eps we get the "Sunny 16" relation between exposure time and ISO;  NULL

t__exp := proc (N, s, R__o) options operator, arrow; H__p(s, theta)/(q*E__ps(N, R__o, theta)) end proc:NULL               evalf(t__exp(16, 100, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)))) = HFloat(0.009683798806264942)*Units:-Unit('s')NULL

t__exp_alt := proc (N, s, R__o) options operator, arrow; k__m*N^2*Pi/(s*q*R__o*E__po) end proc:                  evalf(t__exp_alt(16, 100, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)))) = HFloat(0.00968379880412244)*Units:-Unit('s') 

• 

The Number of Photons NP Reaching the Sensor Area A;

• 

Circle of confusion for 24x36mm "Full Frame" for 1 arcminute view at twice the diagonal:

                          A__cc := Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`^`(Units:-Standard:-`*`(12.6, Unit('`&mu;m`')), 2)), Units:-Standard:-`/`(4)):    

     

• 

  Sensor Bandwidth                                          Photopic Response VP

• 

  Exposure Time for Zone 5: Rscene=12.5% , Saturation in Zone 8 Rscene=100%

• 

  Camera ISO differs from Saturation ISO. Typical Saturation ISO is 2300 when the camera is set to 3200. See DxoMark.

 

NULL

The average number of photons for exposure time based on Reflectance of the scene  relative to the metered value:    

Zone 5;   R__meter := R__scene: 

NP := proc (s, R__o, theta) options operator, arrow; (1/4)*t__exp(N, s, R__meter)*A__cc*q*R__scene*cos(theta)^4*(sum(VP[lambda]*M(lambda*Unit('nm'))*Unit('nm')*lambda*Unit('nm')/(h*c), lambda = 200 .. 4000))/N^2 end proc: 

                                                                               evalf(NP(2300, 1, Units:-Standard:-`*`(0, Unit('deg')))) = HFloat(2191.5645712603696)  NULL

Zone 8;       R__meter := Units:-Standard:-`*`(R__scene, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3))):   NULL

NP__sat := proc (s, theta) options operator, arrow; (1/4)*t__exp(N, s, R__meter)*A__cc*q*R__scene*cos(theta)^4*(sum(VP[lambda]*M(lambda*Unit('nm'))*Unit('nm')*lambda*Unit('nm')/(h*c), lambda = 200 .. 4000))/N^2 end proc:  NULL

                                                                              evalf(NP__sat(2300, Units:-Standard:-`*`(0, Unit('deg')))) = HFloat(17532.516570082957)NULL

NULL

 

Approximate Formula

 

H__sat := proc (s__sat) options operator, arrow; H__p(s__sat, 0)*E__ps(N, 1, 0)/E__ps(N, 1/8, 0) end proc:      

                                                                                       evalf(H__sat(s__sat)) = HFloat(78.53981635)*Units:-Unit(('cd')*('s')/('m')('radius')^2)/s__satNULLNULL

Average Visible Photon Energy

P__e_ave := Units:-Standard:-`*`(Units:-Standard:-`/`(Units:-Standard:-`+`(850, -350)), sum(Units:-Standard:-`*`(Units:-Standard:-`*`(h, c), Units:-Standard:-`/`(Units:-Standard:-`*`(lambda, Unit('nm')))), lambda = 350 .. 850)):                    evalf(P__e_ave) = 0.3533174192e-18*Units:-Unit('J') 

NPtyp := proc (s__sat) options operator, arrow; H__sat(s__sat)*A__cc/(683.002*(Unit('lm')/Unit('W'))*P__e_ave) end proc: 

                               evalf(NPtyp(2300)) = HFloat(17644.363333654386)"(->)"HFloat(17644.363333654386)NULL

NULL

 

Download Photon_Exposure_Array.mw

Obsolete

See my Camera Profiler application instead.

 

This application creates DNG matrices by optimizing Delta E from a raw photo of x-rites color checker. The color temperature for the photograph is also estimated.  Inputs are raw data from RawDigger and generic camera color response from DXO Mark.

Initialization

   

NULL

NULL

NULL

NULL

NULL

XYZoptical to RGB to XYZdata

 

 

Sr,g,b is the relative spectral transmittance of the filter array not selectivity for XY or Z of a given color.

Pulling Sr,g,b out of the integral assumes they are scalars. For example Sr attenuates X, Y and Z by the same amount.

Raw Balance is not White Point Adaptation.

The transmission loss of Red and Blue pixels relative to green is compensated by D=inverse(S). The relation to incident chromaticity, xy is unchanged as S.D=1.

(See Bruce Lindbloom; "Spectrum to XYZ" and "RGB/XYZ Matrices" also, Marcel Patek; "Transformation of RGB Primaries")

 

 

X = (Int(I*xbar*S, lambda))/N:

Y = (Int(I*ybar*S, lambda))/N:

Z = (Int(I*zbar*S, lambda))/N:

N = Int(I*ybar, lambda):

• 

XYZ to RGB

(Vector(3, {(1) = R_Tbb, (2) = G_Tbb, (3) = B_Tbb})) = (Matrix(3, 3, {(1, 1) = XR*Sr, (1, 2) = YR*Sr, (1, 3) = ZR*Sr, (2, 1) = XG*Sg, (2, 2) = YG*Sg, (2, 3) = ZG*Sg, (3, 1) = XB*Sb, (3, 2) = YB*Sb, (3, 3) = ZB*Sb})).(Vector(3, {(1) = X_Tbb, (2) = Y_Tbb, (3) = Z_Tbb}))

NULL

(Vector(3, {(1) = R_Tbb, (2) = G_Tbb, (3) = B_Tbb})) = (Matrix(3, 3, {(1, 1) = Sr, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = Sg, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = Sb})).(Matrix(3, 3, {(1, 1) = XR, (1, 2) = YR, (1, 3) = ZR, (2, 1) = XG, (2, 2) = YG, (2, 3) = ZG, (3, 1) = XB, (3, 2) = YB, (3, 3) = ZB})).(Vector(3, {(1) = X_Tbb, (2) = Y_Tbb, (3) = Z_Tbb}))

 

Camera_Neutral = (Matrix(3, 3, {(1, 1) = Sr, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = Sg, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = Sb})).(Matrix(3, 3, {(1, 1) = XR, (1, 2) = YR, (1, 3) = ZR, (2, 1) = XG, (2, 2) = YG, (2, 3) = ZG, (3, 1) = XB, (3, 2) = YB, (3, 3) = ZB})).(Vector(3, {(1) = X_wht, (2) = Y_wht, (3) = Z_wht}))

NULL

NULL

NULL

• 

RGB to XYZ (The extra step of adaptation to D50 is included below)

 

(Vector(3, {(1) = X_D50, (2) = Y_D50, (3) = Z_D50})) = (Matrix(3, 3, {(1, 1) = XTbbtoXD50, (1, 2) = YTbbtoXD50, (1, 3) = ZTbbtoXD50, (2, 1) = XTbbtoYD50, (2, 2) = YTbbtoYD50, (2, 3) = ZTbbtoYD50, (3, 1) = XTbbtoZD50, (3, 2) = YTbbtoZD50, (3, 3) = ZTbbtoZD50})).(Matrix(3, 3, {(1, 1) = RX*Dr, (1, 2) = GX*Dg, (1, 3) = BX*Db, (2, 1) = RY*Dr, (2, 2) = GY*Dg, (2, 3) = BY*Db, (3, 1) = RZ*Dr, (3, 2) = GZ*Dg, (3, 3) = BZ*Db})).(Vector(3, {(1) = R_Tbb, (2) = G_Tbb, (3) = B_Tbb})) NULL

NULL

(Vector(3, {(1) = X_D50, (2) = Y_D50, (3) = Z_D50})) = (Matrix(3, 3, {(1, 1) = XTbbtoXD50, (1, 2) = YTbbtoXD50, (1, 3) = ZTbbtoXD50, (2, 1) = XTbbtoYD50, (2, 2) = YTbbtoYD50, (2, 3) = ZTbbtoYD50, (3, 1) = XTbbtoZD50, (3, 2) = YTbbtoZD50, (3, 3) = ZTbbtoZD50})).(Matrix(3, 3, {(1, 1) = RX, (1, 2) = GX, (1, 3) = BX, (2, 1) = RY, (2, 2) = GY, (2, 3) = BY, (3, 1) = RZ, (3, 2) = GZ, (3, 3) = BZ})).(Matrix(3, 3, {(1, 1) = Dr, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = Dg, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = Db})).(Vector(3, {(1) = R_Tbb, (2) = G_Tbb, (3) = B_Tbb}))

NULL

(Vector(3, {(1) = X_D50, (2) = Y_D50, (3) = Z_D50})) = (Matrix(3, 3, {(1, 1) = RX_D50, (1, 2) = GX_D50, (1, 3) = BX_D50, (2, 1) = RY_D50, (2, 2) = GY_D50, (2, 3) = BY_D50, (3, 1) = RZ_D50, (3, 2) = GZ_D50, (3, 3) = BZ_D50})).(Matrix(3, 3, {(1, 1) = Dr, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = Dg, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = Db})).(Vector(3, {(1) = R_Tbb, (2) = G_Tbb, (3) = B_Tbb}))

NULL

(Vector(3, {(1) = X_D50wht, (2) = Y_D50wht, (3) = Z_D50wht})) = (Matrix(3, 3, {(1, 1) = RX_D50, (1, 2) = GX_D50, (1, 3) = BX_D50, (2, 1) = RY_D50, (2, 2) = GY_D50, (2, 3) = BY_D50, (3, 1) = RZ_D50, (3, 2) = GZ_D50, (3, 3) = BZ_D50})).(Matrix(3, 3, {(1, 1) = Dr, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = Dg, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = Db})).Camera_Neutral

NULL

Functions

   

NULL

Input Data

   

NULL

Solve for Camera to XYZ D50 and T

   

NULL

 

 

First 46 47 48 49 50 51 52 Last Page 48 of 76