 CoCalc Public FilesTotal Iridium Calculation.sagews
Author: Chris Holden
Views : 38
Compute Environment: Ubuntu 18.04 (Deprecated)
%md
Our goal is to "add up" all of the Iridium deposited on Earth during the major bolide impact 60MYA. We are given a function of the amount of Iridium in ng according to the distance from impact:

$$f(dist)=333.4e^{-0.0008091*dist}.$$

And we look to integrate this function over the surface of the Earth.

The main insight is to think of the function in terms of [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system) ($\rho$, $\theta$, $\phi$)—unfortunately mathematicians and physicists make things more difficult as they have opposite conventions for what symbol represents the vertical angle; here it is $\phi$. Imagine the impact point is the north pole. $\rho$ is constant at 6371, as we are on the surface of the Earth. Distance to the north pole is independent of $\theta$. So distance from that impact point will be determined *only* by the vertical angle $\phi$, starting at 0. Instead of writing

$$f(dist)=333.4e^{-0.0008091*dist},$$

we can write

$$f(\phi)=333.4e^{-0.0008091*6371*\phi}.$$

Something I forgot to write when first documenting my thinking is how $dist$ ends up as $6371*\phi$. It comes down to the designation of radians as units of angular measure: Circumference of an entire circle is $2\pi \rho$ and angular measure around that same span is assigned the value of $2\pi$. So fractions of either angular measure or arc length are computed as the same fraction of $2\pi$.

The job of adding up all of the Iridium across the surface of the Earth is then a surface integral:

$$I = \iint_S{fdS}$$

where $S$ is the surface of the Earth.

It is a standard fact that $dS = \rho^2 sin{\phi}\,d\phi d\theta$. And integrating over the entire Earth means $0 \leq \theta \leq 2\pi$ and $0 \leq \phi \leq \pi$. So the above integral becomes

$$\int_0^{2\pi}\int_0^{\pi}{f(\phi)\rho^2sin{\phi}\,d\phi \,d\theta} = \int_0^{2\pi}\int_0^{\pi}{333.4*e^{-0.0008091*6371*\phi}(6371)^2sin{\phi}\,d\phi \,d\theta}.$$

Note that none of the integrand depends on $\theta$. This makes the integration with respect to $\theta$ trivial, $\int_0^{2\pi}d\theta=2\pi$. And if we let $k=-0.0008091*6371$, then

$$I = 333.4*(6371)^2*2\pi \int_0^{\pi}{e^{k\phi}sin{\phi}\,d\phi}.$$

The integral that is left is not super obvious (I asked the computer to do it before deriving it myself), but it is certainly doable with only a bit of algebra and guessing.

$$\int_0^{\pi}{e^{k\phi}sin{\phi}\,d\phi} = \left[\frac{e^{k\phi}}{k^2+1}(k\sin{\phi}-cos{\phi})\right]_0^{\pi}.$$

Since $\sin{(0)}=\sin{(\pi)}=0$ and $\cos{(0)}=-cos{(\pi)}=1$, this integral evaluates to $\frac{e^{\pi k}+1}{k^2+1}$, and

$$I = 333.4*(6371)^2*2\pi*\frac{e^{-\pi *0.0008091*6371}+1}{(-0.0008091*6371)^2+1}.$$



Our goal is to "add up" all of the Iridium deposited on Earth during the major bolide impact 60MYA. We are given a function of the amount of Iridium in ng according to the distance from impact:

$f(dist)=333.4e^{-0.0008091*dist}.$

And we look to integrate this function over the surface of the Earth.

The main insight is to think of the function in terms of spherical coordinates ($\rho$, $\theta$, $\phi$)—unfortunately mathematicians and physicists make things more difficult as they have opposite conventions for what symbol represents the vertical angle; here it is $\phi$. Imagine the impact point is the north pole. $\rho$ is constant at 6371, as we are on the surface of the Earth. Distance to the north pole is independent of $\theta$. So distance from that impact point will be determined only by the vertical angle $\phi$, starting at 0. Instead of writing

$f(dist)=333.4e^{-0.0008091*dist},$

we can write

$f(\phi)=333.4e^{-0.0008091*6371*\phi}.$

Something I forgot to write when first documenting my thinking is how $dist$ ends up as $6371*\phi$. It comes down to the designation of radians as units of angular measure: Circumference of an entire circle is $2\pi \rho$ and angular measure around that same span is assigned the value of $2\pi$. So fractions of either angular measure or arc length are computed as the same fraction of $2\pi$.

The job of adding up all of the Iridium across the surface of the Earth is then a surface integral:

$I = \iint_S{fdS}$

where $S$ is the surface of the Earth.

It is a standard fact that $dS = \rho^2 sin{\phi}\,d\phi d\theta$. And integrating over the entire Earth means $0 \leq \theta \leq 2\pi$ and $0 \leq \phi \leq \pi$. So the above integral becomes

$\int_0^{2\pi}\int_0^{\pi}{f(\phi)\rho^2sin{\phi}\,d\phi \,d\theta} = \int_0^{2\pi}\int_0^{\pi}{333.4*e^{-0.0008091*6371*\phi}(6371)^2sin{\phi}\,d\phi \,d\theta}.$

Note that none of the integrand depends on $\theta$. This makes the integration with respect to $\theta$ trivial, $\int_0^{2\pi}d\theta=2\pi$. And if we let $k=-0.0008091*6371$, then

$I = 333.4*(6371)^2*2\pi \int_0^{\pi}{e^{k\phi}sin{\phi}\,d\phi}.$

The integral that is left is not super obvious (I asked the computer to do it before deriving it myself), but it is certainly doable with only a bit of algebra and guessing.

$\int_0^{\pi}{e^{k\phi}sin{\phi}\,d\phi} = \left[\frac{e^{k\phi}}{k^2+1}(k\sin{\phi}-cos{\phi})\right]_0^{\pi}.$

Since $\sin{(0)}=\sin{(\pi)}=0$ and $\cos{(0)}=-cos{(\pi)}=1$, this integral evaluates to $\frac{e^{\pi k}+1}{k^2+1}$, and

$I = 333.4*(6371)^2*2\pi*\frac{e^{-\pi *0.0008091*6371}+1}{(-0.0008091*6371)^2+1}.$

We can then use Sage to compute this number approximately.

n(333.4*(6371)^2*2*pi*(exp(-pi*0.0008091*6371)+1)/((-0.0008091*6371)^2+1))

3.08387600430190e9

Unfortunately, this isn't the number we want. We were hoping for something closer to $1\times 10^{19}$ ng. Certainly I could have made many errors. Most of the computations here, the number crunching and the finding of tricky integrals, I did in two separate ways and got the same answer, but I'm not the most historically accurate data entry clerk. I tried to spell everything out as clearly as possible so mistakes would be easier to find.

The earlier error, where we had a negative value for $k$ does not affect the problem much analytically. In general, we should expect a smaller answer than with $k$ positive, as the resulting unsigned exponent is larger than one. And in fact that is what we got. Instead of $3\times 10^{16}$ we now get $3\times 10^9$.

Just to check that I haven't misplaced a parenthesis or something similar, let's do the calculation in parts.

First, the constant out in front:

C = n(333.4*(6371)^2*2*pi)
C

8.50277474673617e10

Now the exponent

E = n(exp(-pi*0.0008091*6371*10^6))
E

2.39778883022791e-7033055

And now $k^2+1$:

D = n((-0.0008091*6371)^2+1)
D

27.5717166411312

Putting it all together

C*(E+1)/D

3.08387571851504e9