CoCalc Public FilesTotal Iridium Calculation.sagewsOpen with one click!
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.4e0.0008091dist.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.4e0.0008091dist,f(dist)=333.4e^{-0.0008091*dist},

we can write

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

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

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

I=SfdSI = \iint_S{fdS}

where SS is the surface of the Earth.

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

02π0πf(ϕ)ρ2sinϕdϕdθ=02π0π333.4e0.00080916371ϕ(6371)2sinϕdϕdθ.\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, 02πdθ=2π\int_0^{2\pi}d\theta=2\pi. And if we let k=0.00080916371k=-0.0008091*6371, then

I=333.4(6371)22π0πekϕsinϕdϕ.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.

0πekϕsinϕdϕ=[ekϕk2+1(ksinϕcosϕ)]0π.\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(π)=0\sin{(0)}=\sin{(\pi)}=0 and cos(0)=cos(π)=1\cos{(0)}=-cos{(\pi)}=1, this integral evaluates to eπk+1k2+1\frac{e^{\pi k}+1}{k^2+1}, and

I=333.4(6371)22πeπ0.00080916371+1(0.00080916371)2+1.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×10191\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 kk does not affect the problem much analytically. In general, we should expect a smaller answer than with kk positive, as the resulting unsigned exponent is larger than one. And in fact that is what we got. Instead of 3×10163\times 10^{16} we now get 3×1093\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 k2+1k^2+1:

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

Putting it all together

C*(E+1)/D
3.08387571851504e9