## Nicole Sharp

Miss Nicole Sharp

## 130 Reputation

1 years, 150 days
Frostburg State University (FSU)
Kappa Mu Epsilon Alumna
Cumberland, Maryland, United States

## Social Networks and Content at Maplesoft.com

Nicole Sharp of Frostburg State University (FSU) and Allegany College of Maryland (ACM), United States of America (USA). https://www.nicolesharp.net/

## @dharr, not sure if you're using the...

@dharr, not sure if you're using the same eccentricity.  The eccentricity I got doesn't work:

```AddConstant(Varda_diameter, symbol = d[Varda], value = 740., uncertainty = 14., units = km) :

AddConstant(Varda_radius, symbol = r[Varda], derive = d[Varda]/2) :

# f[Varda] := Quantity(0.080, 0.049) :
AddConstant(Varda_flattening, symbol = f[Varda], value = 0.080, uncertainty = 0.049, units = 1) :

AddConstant(Varda_eccentricity, symbol = e[Varda], derive = sqrt(2*f[Varda] - f[Varda]^2)) :

AddConstant(Varda_equatorial_radius, symbol = r[e,Varda], derive = 3*r[Varda]*sqrt(1 - e[Varda]^2)) :
AddConstant(Varda_polar_radius, symbol = r[p,Varda], derive = 3*r[Varda]*sqrt(1 - e[Varda]^2)/(2 + sqrt(1 - e[Varda]^2))) :

# AddConstant(Varda_equatorial_radius, symbol = r[e,Varda], derive = -3*r[Varda]/(f[Varda] - 3)) :
# AddConstant(Varda_polar_radius, symbol = r[p,Varda], derive = 3*(f[Varda] - 1)*r[Varda]/(f[Varda] - 3)) :```

Using this has the advantage though that both the equatorial and polar radii are now dimensionless instead of only the polar radius being dimensionless.  This gives a dimensionless quantity for both area and volume, but density is in kilograms instead of kilograms per cubic meter.

Unfortunately this still creates serious problems.  For example, you can't define a derivative constant with an average density if one of the densities is in kilograms instead of kilograms per cubic meter.  Will just have to not use the ScientificConstants package until this can get fixed.

I also noticed that "with(Statistics) : Mean(%)" doesn't seem to work with Units either.

```d[Varda] := Quantity(740., 14.)*Unit(km) :
f[Varda] := Quantity(0.080, 0.049) :
m[Varda] := Quantity(2.45E20, 0.06E20)*Unit(kg) :

r[Varda] := d[Varda]/2 :

r[e,Varda] := -3*r[Varda]/(f[Varda] - 3) :
r[p,Varda] := 3*(f[Varda] - 1)*r[Varda]/(f[Varda] - 3) :

V[Varda] := V_spheroid(r[e,Varda], r[p,Varda]) :

rho[Varda] := m[Varda]/V[Varda] :

rho[cubewano] := (rho[Varda] + Constant(rho[Varuna], units))/2 :

# AddConstant(cubewano_density, symbol = rho[cubewano], derive = (rho[Varda] + rho[Varuna])/2) :
```

A better example of how serious this bug is:

```# f[Varda] := Quantity(0.080, 0.049) :
AddConstant(Varda_flattening, symbol = f[Varda], value = 0.080, uncertainty = 0.049) :

AddConstant(Varda_diameter, symbol = d[Varda], value = 740., uncertainty = 14., units = km) :

AddConstant(Varda_radius, symbol = r[Varda], derive = d[Varda]/2) :

AddConstant(Varda_equatorial_radius, symbol = r[e,Varda], derive = -3*r[Varda]/(f[Varda] - 3)) :
AddConstant(Varda_polar_radius, symbol = r[p,Varda], derive = 3*(f[Varda] - 1)*r[Varda]/(f[Varda] - 3)) :
```

"AddConstant" will not accept "Quantity" for "derive".  So there is no way to calculate the area or volume of this spheroid with error and units unless not using the "ScientificConstants" package, not using units, or using the eccentricity workaround.

What is most strange is that the polar radius is dimensionless but the equatorial radius has the correct units (meters).  Unfortunately this means that the volume is in square meters instead of cubic meters when multiplying the dimensioned length by the dimensionless length.

## @dharr, unfortunately that doesn't w...

Unfortunately that doesn't seem to work either.  I opened a new thread for this new issue with integrating on Maple:

https://www.mapleprimes.com/questions/237724-Ellipsoids-With-Units-And-Error?sq=237724

## @Nicole Sharp, a similar problem with no...

A similar problem with non-cancelable units appears to occur when using the "ellipsoid" function (which uses elliptic integrals) together with Units and ScientificErrorAnalysis.

```alias(asin = arcsin) :
alias(pi = Pi) :
alias(S_ellipsoid = ellipsoid) :
E := (n) -> 10^n :
with(ScientificErrorAnalysis) :
with(Units) :

_km := Unit(km) :
_lm := Unit(lm) :
_lx := Unit(lx) :
_m := Unit(m) :
_m2 := Unit(m^2) :
_sr := Unit(sr) :
AddUnit(nit, symbol = nt, prefix = SI) :
f_Terra := 1/Quantity(298.257222101, 0.000000001) :
R_Io := Quantity(0.63, 0.02) :
R_Terra := Quantity(0.306, 0.001) :
S_ellipse := (r_e, r_p) -> pi*r_e*r_p :
S_spheroid := (r_e, r_p) -> S_ellipsoid(r_e, r_e, r_p) :

_nt := Unit(nt) :
a_Terra := Quantity(149598023., 1.)*_km :
alpha_Jupiter := Quantity(816.363, 0.001)*E(6)*_km :
A_Terra := 1 - R_Terra :
d_e_Io := Quantity(3660.0, 0.1)*_km :
d_m_Io := Quantity(3637.4, 0.1)*_km :
pi_Jupiter := Quantity(740.595, 0.001)*E(6)*_km :
Phi_Sol := Quantity(3.75, 0.01)*E(28)*_lm :
r_e_Terra := Quantity(6378.1370, 0.0001)*_km :
r_mu_Io := Quantity(1821.6, 0.5)*_km :
S_sphere := (r) -> S_spheroid(r, r) :

a_Jupiter := (alpha_Jupiter + pi_Jupiter)/2 :
r_e_Io := d_e_Io/2 :
r_m_Io := d_m_Io/2 :
r_p_Terra := -(f_Terra - 1)*r_e_Terra :

E_Sol_Io := Phi_Sol/S_sphere(a_Jupiter) :
r_e_mu_Io := (r_e_Io + r_m_Io)/2 :
r_mu_Terra := (2*r_e_Terra + r_p_Terra)/3 :
r_p_Io := 3*r_mu_Io - r_e_Io - r_m_Io :

Delta_Terra := r_mu_Terra :
M_Io := R_Io*E_Sol_Io :
S_Io := S_ellipsoid(r_e_Io, r_m_Io, r_p_Io) :

Delta_Jupiter := sqrt(a_Jupiter^2 + a_Terra^2) - Delta_Terra :
Phi_Io := M_Io*S_Io/2 :

Delta_Io := Delta_Jupiter :

E_Io := A_Terra*Phi_Io/S_sphere(Delta_Io) :
rho_e_mu_Io := asin(r_e_mu_Io/Delta_Io)*_rad :
rho_p_Io := asin(r_p_Io/Delta_Io)*_rad :

Omega_Io := S_ellipse(rho_e_mu_Io, rho_p_Io) :

L_Io := E_Io/Omega_Io : # Ionian luminance in nits (candelas per square meter)

combine(combine(L_Io, units), errors);```

## I think it might be a bug with the Units...

The Units package seems buggy.  As far as I can tell, it hasn't been updated since at least 2011, since the definition of the astronomical unit (AU) is incorrect.  There doesn't seem to be a way to remove units from this dimensionless quantity below.  Using "convert(%, unit_free)" doesn't work either.  Nor does "convert(%, units, 1)".

```alias(Gamma = GAMMA) :
alias(pi = Pi) :
alias(zeta = Zeta) :
E := (n) -> 10^n :
with(ScientificErrorAnalysis) :
with(Units) :

_Hz := Unit(Hz) :
_J := Unit(J) :
_K := Unit(K) :
_m := Unit(m) :
_m2 := Unit(m^2) :
_nm := Unit(nm) :
_s := Unit(s) :
_W := Unit(W) :

c := 299792458*_m/_s :
h := 662607015*E(-8)*E(-34)*_J/_Hz :
k := 1380649*E(-6)*E(-23)*_J/_K :
lambda_L_max := Quantity(750., 10.)*_nm :
lambda_L_min := Quantity(380., 10.)*_nm :
T_e_Sol := Quantity(5772.0, 0.8)*_K :

M_L_Sol := pi*Int(2*c^2*h/((exp(c*h/(k*lambda*T_e_Sol)) - 1)*lambda^5), lambda = lambda_L_min .. lambda_L_max, method = _Dexp) :
sigma := 2*pi*Gamma(4)*zeta(4)*k^4/(c^2*h^3) :

M_Sol := T_e_Sol^4*sigma :

eta_L_Sol := M_L_Sol/M_Sol :

combine(combine(eta_L_Sol, errors), units);```

## @dharr, I agree that this seems like a b...

@dharr, this seems like a bug when it gives a different answer for the same integration depending on whether or not units are included.  Everything is in SI units so there are no unit conversions (all of the units used equal one) so it should have exactly the same result both with and without units.

"method = _Dexp" is not required when doing this integration on Maple for Excel when the Excel values are saved without units and the Excel precision is limited to 15 significant figures.

## @C_R, thank you!!I went through all of t...

@C_R, thank you!!

I went through all of the Maple 2023 Help pages I could find that might be relevant and I actually did try "method = _Dexp" yesterday but it gave me an error message.  I didn't realize that I also need to change "int" to "Int".  I just tried it again and it works now!

Whew.  That took a few days of head-pounding but I am glad it appears to be working correctly now.

I still think that this should be considered as a bug with Maple.  There should be an option in the Maple settings to instruct Maple that all integrations should only give real results and thus Maple can then automatically try alternative integration methods whenever a non-real result is detected.  This appears to be what Wolfram Alpha does since Mathematica does not give a complex number for this integral.

```combine(Pi*Int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = lambda_V_min .. lambda_V_max, method = _Dexp), errors);

# Quantity(2.7551199571700602410253741952575E7, 17871.934324613671637552924319647)*_W/_m2```

## I was able to get the error from the int...

I was able to get the error from the integral but only as a complex number.  "fnormal" doesn't appear to work here and attempting "evalf" will remove the error.  Any suggestions?

The imaginary part is no larger than E-22*I but even setting "fnormal(%, 4)" doesn't work.

```alias(Gamma = GAMMA) :
alias(pi = Pi) :
alias(zeta = Zeta) :
E := (n) -> 10^n :
with(ScientificErrorAnalysis) :
with(Units) :
with(ScientificConstants) :

_Hz := Unit(Hz) :
_J := Unit(J) :
_K := Unit(K) :
_m := Unit(m) :
_m2 := Unit(m^2) :
_nm := Unit(nm) :
_s := Unit(s) :
_W := Unit(W) :
AddConstant(Solar_temperature, symbol = T[Sol], value = 5772.0, uncertainty = 0.8, units = K) :

c := 299792458*_m/_s :
h := 662607015*E(-8)*E(-34)*_J/_Hz :
k := 1380649*E(-6)*E(-23)*_J/_K :
lambda_V_max := 750*_nm :
lambda_V_min := 380*_nm :
T_Sol := Constant(T[Sol], units) :
T0_Sol := 5772*_K :

M := (T, lambda_min, lambda_max) -> pi*int(2*c^2*h/((exp(c*h/(k*lambda*T)) - 1)*lambda^5), lambda = lambda_min .. lambda_max) :
sigma := 2*Gamma(4)*pi*zeta(4)*k^4/(c^2*h^3) :

M_inf := (T) -> T^4*sigma :

M_Sol := M_inf(T_Sol) :
M0_Sol := M_inf(T0_Sol) :

(*

combine(M_inf(Constant(T[Sol], units)), errors);
# Quantity(6.2938592470335950467548474587302E7, 34893.190558744809684018558329758)*_W/_m2

simplify(fnormal(combine(pi*int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = lambda_V_min .. lambda_V_max), errors)), zero);
# Quantity(2.7551199571700602410253741952750E7 - 6.9115038378975451246178154432149E-22*I, 17871.934324613671637552924319095 + 2.2619467105846511316931032359613E-25*I)*_W/_m2

*)```

## Still trying to find something that work...

Still trying to find something that works consistently but the integral seems more well-behaved when the minimum and maximum wavelengths are not exact.

Best I have so far is:

`simplify(fnormal((Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = evalf(lambda_V_min) .. evalf(lambda_V_max))),16),zero);`

The Maple 2023 Help does provide a function to give different integration settings depending on whether inputting exact or inexact values, outlined below.

"You can define a new function and corresponding evalf subroutine such that the evalf subroutine is called if the arguments are of floating-point type. To clearly show which code Maple uses, in the following example, the mathematical function returns 2 times the argument for non-floating-point input and 2.01 times the argument for floating-point input."

NewFunction := proc(x)
local s;
if type(x,'complex'('float')) then
s := evalf('NewFunction'(x));
if type(s,'complex'('float')) then
return s;
end if;
end if;
return 2*x;
end proc:
`evalf/NewFunction` := proc(xx)
local x;
x := evalf(xx);
if not type(x, 'complex'('float')) then
return 'NewFunction'(x);
end if;
return 2.01*x;
end proc:

## Here is something that works with units ...

Here is something that works with units but the error is still removed.

```simplify(fnormal(combine(Pi*int(2*(299792458*Unit(m)/Unit(s))^2*662607015*10^(-8)*10^(-34)*(Unit(J)/Unit(Hz))/((exp(299792458*(Unit(m)/Unit(s))*662607015*10^(-8)*10^(-34)*(Unit(J)/Unit(Hz))/(1380649*10^(-6)*10^(-23)*(Unit(J)/Unit(K))*lambda*Quantity(5772.0, 0.8)*Unit(K))) - 1)*lambda^5), lambda = 380*10^(-9)*Unit(m) .. 750*10^(-9)*Unit(m), numeric = true), errors)), zero);
```

"simplify(fnormal(combine(%,errors)),zero)" should provide "Quantity(x, y)*Unit(z)".

## @dharr, the behavior in Maple is still v...

@dharr, the behavior in Maple is still very inconsistent.  Here are four examples.

`Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772.0)) - 1)*lambda^5), lambda = 380*10^(-9) .. 750*10^(-9));`

Using "5772.0" as a dimensionless quantity-without-error produces the correct real result.

`Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*Quantity(5772.0,0.8))) - 1)*lambda^5), lambda = 380*10^(-9) .. 750*10^(-9));`

Using "Quantity(5772.0, 0.8)" as a dimensionless quantity-with-error produces the correct real result but removes the error.  "combine(%, errors)" does not work.

`Pi*int(2*(299792458*Unit(m)/Unit(s))^2*662607015*10^(-8)*10^(-34)*(Unit(J)/Unit(Hz))/((exp(299792458*(Unit(m)/Unit(s))*662607015*10^(-8)*10^(-34)*(Unit(J)/Unit(Hz))/(1380649*10^(-6)*10^(-23)*(Unit(J)/Unit(K))*lambda*5772.0*Unit(K))) - 1)*lambda^5), lambda = 380*10^(-9)*Unit(m) .. 750*10^(-9)*Unit(m));`

Using "5772.0*Unit(K)" produces a complex number but with the correct units.

`Pi*int(2*(299792458*Unit(m)/Unit(s))^2*662607015*10^(-8)*10^(-34)*(Unit(J)/Unit(Hz))/((exp(299792458*(Unit(m)/Unit(s))*662607015*10^(-8)*10^(-34)*(Unit(J)/Unit(Hz))/(1380649*10^(-6)*10^(-23)*(Unit(J)/Unit(K))*lambda*Quantity(5772.0,0.8)*Unit(K))) - 1)*lambda^5), lambda = 380*10^(-9)*Unit(m) .. 750*10^(-9)*Unit(m));`

Using "Quantity(5772.0, 0.8)*Unit(K)" produces a complex number with the correct units but the error is again removed from the answer.

So there appear to be two problems here.  Somehow adding units causes a complex number to be returned but keeping the quantities dimensionless returns a real number instead.  Then there is another problem that the error appears to be removed.  I don't know if ScientificErrorAnalysis is designed to work with these types of integrations but there should be a way to calculate the error since there is only one input with error (all other values used are exact quantities).

I will try fiddling with "fnormal" and "simplify" to see if I can get something that might work consistently.  Otherwise another option might be to remove the units before integrating and then add them back after integration.  It is frustrating though since this does not seem to be a problem on Wolfram Alpha.  According to Maplesoft, Maple is supposed to be better than Mathematica at this kind of thing.

## with(Units) :c := 299792458*Unit(m...

"5772.0*Unit(K)" must be integrated symbolically for wavelengths from 0 to infinity.

"5772.0*Unit(K)" must be integrated numerically for wavelengths from 380 nm to 750 nm.

"5772*Unit(K)" can be integrated either symbolically or numerically for wavelengths from 0 to infinity.

"5772*Unit(K)" must be integrated numerically for wavelengths from 380 nm to 750 nm.

I will have to experiment around, but the numeric integration might be the best option for finite bandwidths and then the symbolic integration might only be required for infinite bandwidth (which is equal to the algebraic Stefan-Boltzmann law without needing integration).

```with(Units) :

c := 299792458*Unit(m)/Unit(s) :
h := 662607015*E(-8)*E(-34)*Unit(J)/Unit(Hz) :
k := 1380649*E(-6)*E(-23)*Unit(J)/Unit(K) :
lambda_V_max := 750*Unit(nm) :
lambda_V_min := 380*Unit(nm) :
T_Sol := 5772.0*Unit(K) :
T0_Sol := 5772*Unit(K) :

(*

# radiant exitances for non-exact temperatures:
# evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = 0*_m .. infinity*_m, numeric = true));
evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = 0*_m .. infinity*_m, numeric = false)); # 6.2938592470335950467548474587301E7*_W/_m2
evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = lambda_V_min .. lambda_V_max, numeric = true)); # 2.7551199571700602410253741952570E7*_W/_m2
# evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = lambda_V_min .. lambda_V_max, numeric = false));

# radiant exitances for exact temperatures:
evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T0_Sol)) - 1)*lambda^5), lambda = 0*_m .. infinity*_m, numeric = true)); # 6.2938592470335950467548474587304E7*_W/_m2
evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T0_Sol)) - 1)*lambda^5), lambda = 0*_m .. infinity*_m, numeric = false)); # 6.2938592470335950467548474587300E7*_W/_m2
evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T0_Sol)) - 1)*lambda^5), lambda = lambda_V_min .. lambda_V_max, numeric = true)); # 2.7551199571700602410253741952570E7*_W/_m2
# evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T0_Sol)) - 1)*lambda^5), lambda = lambda_V_min .. lambda_V_max, numeric = false));
*)```

## @dharr, numeric integration works from z...

@dharr, numeric integration works from zero to infinity for "5772*Unit(K)" but does not work for "5772.0*Unit(K)".  If using a non-exact value, then it requires symbolic integration instead.  This is part of the unusual behavior I don't like.

It would be nice if there was an option for "evalf" where you could say:

"evalf(int(%), real)"

where:

if evalf(int(5772)) is complex

then evalf(int(5772, numeric = true))

and if evalf(int(5772.0)) is complex

then evalf(int(5772.0, numeric = false))

That should automatically always give the real result, like Wolfram Alpha does.

Otherwise it is impossible to create a generic function, since the integration has to be changed depending on the inputted values.

## numeric or symbolic...

@dharr, using "numeric" seems like the simplest solution.  Except that there isn't an intuitive way to know when to use this, as shown in the examples below.  The integration from 380 nm to 750 nm requires numeric integration for both "5722 K" and "5772.0 K" whereas the integration from 0 to infinity requires symbolic integration for both "5772 K" and "5772.0 K".

What is interesting is that Wolfram Alpha does not have this problem.  Wolfram Alpha always returns a real answer, not a complex one, without the user needing to specify whether to use symbolic or numeric integration.  Is there a way to instruct Maple to automatically choose either symbolic or numeric integration, based on whichever one will return a real (noncomplex) result for "evalf"?  Otherwise the user is forced through trial and error on each new integration to find which method gives real answers.

```evalf(Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772.0)) - 1)*lambda^5), lambda = 0 .. infinity, numeric = false));
7
6.2938592470335950467548474587301 10

evalf(Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772)) - 1)*lambda^5), lambda = 0 .. infinity));
7
6.2938592470335950467548474587300 10

evalf(Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772.0)) - 1)*lambda^5), lambda = 380*10^(-9) .. 750*10^(-9)));
7
2.7551199571700602410253741952570 10

evalf(Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772)) - 1)*lambda^5), lambda = 380*10^(-9) .. 750*10^(-9), numeric));
7
2.7551199571700602410253741952570 10 ```

## 5772.0...

Seems to be the opposite problem here.  In this case, using "5772" returns a complex result but using "5772.0" returns a real result.  I still don't like this type of behavior on Maple though, especially when Wolfram Alpha always returns the correct real result.  Is there another setting in Maple to fix this?

 1 2 3 4 Page 1 of 4
﻿