In 1938, physicist Frank Benford made a fascinating discovery that would later revolutionize fraud detection and data analysis. He noticed that in naturally occurring numbers, the first digit isn't randomly distributed as one might expect. Instead, it follows a peculiar logarithmic pattern now known as Benford's Law.

**In 1938, physicist Frank Benford made a fascinating discovery that would later revolutionize fraud detection and data analysis. He noticed that in naturally occurring numbers, the first digit isn't randomly distributed as one might expect. Instead, it follows a peculiar logarithmic pattern now known as Benford's Law.**

## Benford's Law Formula

According to this law, in many real-world datasets, the number 1 appears as the first digit about 30.1% of the time, while 9 appears only 4.6% of the time. This counter-intuitive distribution shows that lower digits are far more likely to appear as the first digits than higher ones. The complete distribution predicts that digits 1 through 9 appear as first digits with the following frequencies: 30.1%, 17.6%, 12.5%, 9.7%, 7.9%, 6.7%, 5.8%, 5.1%, and 4.6% according to the following formula:

$$

P(d) = \log_{10}\left( 1 + \frac{1}{d} \right),

$$

where $P(d)$ is the probability of leading digit $d$ to occur.

## Benford's Law Explanation

But why does this pattern emerge? The explanation lies in the way natural processes grow. Consider a city's population that grows by a fixed percentage each year. To increase from 1,000 to 2,000 residents (starting with digit 1 to starting with digit 2), the population needs to double. However, to go from 9,000 to 10,000 residents, it only needs to increase by about 11%. This means numbers tend to spend more time starting with lower digits during natural growth processes.

This mathematical principle appears in an astonishing variety of datasets: stock prices, river lengths, building heights, molecular weights, and even social media statistics. However, Benford's Law only applies to natural, unmanipulated data spanning several orders of magnitude. Data that is artificially constrained or manipulated often deviates from this pattern.

This property makes Benford's Law a powerful tool for detecting potential fraud or data manipulation. When people fabricate numbers, they tend to distribute digits more uniformly, unconsciously creating patterns that deviate from Benford's Law. This principle has been successfully used in forensic accounting, election auditing, and scientific fraud detection.

## Benford's Law Example

Let's see Benford's Law in action using a real-world example. Let us analyze Adobe's [10-K](https://www.sec.gov/ix?doc=/Archives/edgar/data/0000796343/000079634325000004/adbe-20241129.htm) for the fiscal year 2024 using a Chrome extension called [Benford's Law Analyzer](https://chromewebstore.google.com/detail/benfords-law-analyzer/blgcepkhpimpjnanokccgmadcimfhdlj). The extension finds all number sequences in the financial report and automatically calculated their first-digit distribution:

![Benford's law for Adobe 10-K no filtering](https://i.imgur.com/4kbpRNA.png "Benford's law for Adobe 10-K no filtering")

This does not seem to follow Benford's law whatsoever and having in mind what is mentioned above this could be setting off alarm bells. Is Adobe committing financial fraud right before our eyes? Well if we take a closer look at the 10-K we see that a big number of the number sequences is years like 2023 and 2024. In the Chrome extension, we can filter out number sequences and if we filter out all years from 2016 to 2026 we instead get this distribution:

![Benford's law for Adobe 10-K no filtering](https://i.imgur.com/jVEkS72.png "Benford's law for Adobe 10-K no filtering")

Now the distribution is looking more like what you would expect although it is not perfect, but the reason for this might be another logical explanation like that of the years. The lesson here is that you need to know and understand the dataset you are working with and if there are number sequences that might be occurring at an unusual rate compared to the rest of the dataset.

This simple yet powerful principle continues to find new applications in our data-driven world. From detecting financial fraud to validating scientific research, Benford's Law is a quick and effective way to spot potential anomalies in numerical data. While deviation from Benford's Law doesn't necessarily prove manipulation - there are many legitimate reasons why data might not follow this pattern - it provides a valuable first step in identifying datasets that might deserve closer scrutiny.

Benford's Law - Analysis and Example

Linear and angular momentum has certain properties. Here we explore momentum's equations, rate of change, conservation, and link to torque and impulses.

**Linear and angular momentum has certain properties. Here we explore momentum's equations, rate of change, conservation, and link to torque and impulses.**

## Momentum Equation

The formula for momentum is

$$
    p = mv,
$$

where $m$ is the mass of a moving object with velocity $v$. In vector form it becomes

$$
    \vec{p} = m \vec{v}.
$$

## Change in Momentum Formula

To calculate the change in momentum you have to use the change in velocity $\Delta \vec{v} = \vec{v}_2 - \vec{v}_1$ and plug it into the expression to get the change in momentum:

$$
    \Delta \vec{p} = m \Delta \vec{v}
$$

## Rate of Change of Momentum

To get the rate of change of momentum we will take the time derivative of the momentum equation:

$$
    \frac{d \vec{p}}{dt} = \frac{d}{dt} (m \vec{v}).
$$

Since the mass is assumed to be constant the derivative only acts on the velocity and the derivative of velocity is acceleration ($\vec{a}$) so we get:

$$
    \frac{d \vec{p}}{dt} = m \vec{a} = \Sigma \vec{F}.
$$

As you can see we get the expression for Newton's second law which says that the sum of all forces acting upon an object is equal to the mass times acceleration of the body.

## Is Impulse The Change in Momentum?

To get the impulse we have to integrate the equation above over time. In other words, we have to integrate the forces acting on an object over time:

$$
    \int_{t_1}^{t_2} \Sigma \vec{F} dt = \int_{t_1}^{t_2} \frac{d \vec{p}}{dt} dt = \vec{p}_2 - \vec{p}_1 = \Delta \vec{p},
$$

where $\vec{p}_1$ and $\vec{p}_2$ is the linear momentum at time $t_1$ and $t_2$ respectively. As can be seen, the impulse is the change in momentum.

## Conservation of Momentum Equation

Observing the equation above, if the resultant force $\Sigma \vec{F}$ is zero then the momentum does not change and is conserved. Now consider two objects, $A$ and $B$, colliding. We know from Newton's third law that each action (force) has an equal opposite reaction. So the change in momentum for object $A$ is equal to the negative of the change of object $B$, $\Delta \vec{p}_1 = -\Delta \vec{p}_2$. So the total change for the system is

$$
    \Delta \vec{p} = \Delta \vec{p}_1 + \Delta \vec{p}_2 = 0
$$

and the momentum is conserved. 

## Angular Momentum Equation

The formula for angular momentum $\vec{L}_O$ for an object with mass $m$ around origin $O$ is

$$
    \vec{L}_O = \vec{r} \times m \vec{v},
$$

where $\vec{r}$ is the distance in vector form between the origin and the object and $\vec{v}$ is the object's velocity. If it is broken down into components then you get

$$
    L_x = m (v_z y - v_y z)
$$
$$
    L_y = m (v_x z - v_z x)
$$
$$
    L_z = m (v_y x - v_x y)
$$

## Angular Momentum Units

Looking at the equation above the angular momentum has the units $\text{kg} \cdot \text{m}^2/\text{s} = \text{N} \cdot \text{m} \cdot \text{s}$ and in US customary units that is lb-ft-sec.

## Rate of Change of Angular Momentum

To get the rate of change of the angular momentum we have to take the time derivative of the equation above which gives us

$$
    \frac{d \vec{L}_O}{dt} = \frac{d \vec{r}}{dt} \times m \vec{v} + \vec{r} \times m \frac{d \vec{v}}{dt}
$$
$$
    = \vec{v} \times m \vec{v} + \vec{r} \times m \vec{a}.
$$

The cross product of parallel vectors is zero so $\vec{v} \times m \vec{v}$ is zero. Therefore we have the following equation for the angular momentum rate of change:

$$
    \frac{d \vec{L}_O}{dt} = \vec{r} \times m \vec{a}.
$$

For the linear momentum, $\Sigma \vec{F}$ represented the resultant of all forces acting upon the object, and with Newton's law, we know that $\Sigma \vec{F} = m \vec{a}$. So for a total moment $\Sigma \vec{M}_O$ around the origin $O$ is given by

$$
    \Sigma \vec{M}_O = \vec{r} \times m \vec{a}.
$$

As you can see the formula is the same as the one for the rate of change. In other words, the rate of change of angular momentum is equal to the total moment:

$$
    \Sigma \vec{M}_O = \frac{d \vec{L}_O}{dt}.
$$

## Relationship Between Torque And Angular Momentum

Since torque is a moment that causes rotation around an axis the formula above also gives us the formula for the torque around a fixed point $O$:

$$
    \vec{\tau} = \frac{d \vec{L}_O}{dt}.
$$

## Angular Impulse And Angular Momentum

Now that we know that the total moment acting upon an object is equal to the rate of change of the angular momentum we can obtain the effect of the moment on the angular momentum over a time interval using this formula:

$$
    \int_{t_1}^{t_2} \Sigma \vec{M}_O dt = \int_{t_1}^{t_2} \frac{d \vec{L}_O}{dt} dt
$$
$$
    = [\vec{L}_O]_2 - [\vec{L}_O]_1 = \Delta \vec{L}_O.
$$

The angular impulse is defined as the product of moment and time, so the formula above states that the total angular impulse is equal to the change in angular momentum.

## Conservation of Angular Momentum Equation

Similar to the conservation of linear momentum if the resultant moment about a fixed point $O$ is zero over an interval of time then angular momentum $\vec{L}_O$ does not change and is conserved, as seen in the equation above. If two objects, $A$ and $B$, interact over an interval of time then we know from Newton's third law that each action (force) has an equal opposite reaction and this also applies to the moment, so $\Delta \vec{L}_A = -\Delta \vec{L}_B$. The total change in the angular momentum of the system is

$$
    \Delta \vec{L} = \Delta \vec{L}_A + \Delta \vec{L}_B = 0
$$

and we can see that the total angular momentum is conserved.

Momentum Equation - Linear And Angular Conservation

The continuity equation for fluids can be expressed in cartesian, cylindrical, or spherical coordinates, each with distinct properties and formulas.

**The continuity equation for fluids can be applied to problems like calculating the outlet velocity of a pipe with varying diameters. It can also be expressed in different coordinate systems such as the cylindrical and spherical coordinate systems. Let us explore the properties of the different formulas for the continuity equation for fluids.**

## Continuity Equation For Steady Flow

Let us assume a fluid moves in a steady flow through a pipe as in the figure above, meaning the mass flow ($\dot{m}$) into $A_1$ is the same as the mass flow out of $A_2$. This leads to the following formula:

$$
    \dot{m} = \rho_1 v_1 A_1 = \rho_2 v_2 A_2,
$$

where $v_1$ and $v_2$ are the inlet and outlet velocities respectively. $A_1$ is the area of the inlet and $A_2$ is the area of the outlet. $\rho_1$ is the density of the fluid at the inlet and $\rho_2$ is the density of fluid at the outlet. To understand the formula you can think of the velocity $v$ as a height per second and $A$ as a base area so when the two multiply you get a volume per second. $\rho$ is the density so that gives us mass per second, aka the mass flow.

### Continuity Equation Compressible Flow

If the fluid is compressible meaning that the density changes when the pressure changes then the equation above should be used since when the area, as an example, decreases the density and/or velocity have to increase to maintain the equality in the formula.

### Continuity Equation Incompressible Flow

If the fluid is incompressible, which water essentially is, then density is constant ($\rho_1 = \rho_2$) so we can cancel it out on both sides and get:

$$
    v_1 A_1 = v_2 A_2.
$$

So if the area decreases the velocity has to increase. You most likely have experienced this effect when putting your thumb over the end of a water hose. When you do that the area decreases and you get a high velocity flow exiting the hose.

## Continuity Equation in Differential Form

The continuity equation in differential form for fluids looks like this:

$$
    \frac{\partial \rho}{\partial t} = - \nabla \cdot (\rho \vec{u}),
$$

where $\rho$ is the density, $t$ is time, and $\vec{u}$ the flow velocity vector field. If you would like a derivation of this formula then check out our post on the [continuity equation in electromagnetism](https://www.engineeringguidance.com/continuity-equation-in-electromagnetism) and just replace the charge density $\rho$ with the fluid density and the current density $\vec{j}$ with $\rho \vec{u}$.

## Continuity Equation in Cartesian Coordinates

The continuity equation in cartesian coordinates $(x,y,z)$ is

$$
    \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} + \frac{\partial (\rho w)}{\partial z} = 0,
$$

where $u$, $v$, and $w$ are the flow velocity in the $x$-, $y$-, and $z$-direction respectively. 

### Continuity Equation in Cartesian Coordinates With Incompressible Flow

If the fluid is incompressible, $\rho$ is constant then the formula becomes:

$$
    \rho \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} \right) = 0,
$$

since $\frac{\partial \rho}{\partial t}=0$ (<a href="https://www.amazon.com/Handbook-Mathematics-I-N-Bronshtein/dp/3662462206?crid=2EY5M5H7FT2IH&dib=eyJ2IjoiMSJ9.DwfzsUXXkM7e94_ArqNkD7rI9LmiHe7Nl2zScHPEeCKOU1geuDYJusAIGYH94tYi0OQ9FNiQ4yi64-rEENwjhruC-ruHygf1MlPUoylnNrhbsTbUckL03fNH60wBfKaC-Uwew6hftZIJpPJZUSt1u1-MP39iufzwGsdI5WdbGMJm9l7J_Ea3MFirPBITTnUrW7Zx80gGJl0GXL1JdcZzsvBPWQ944WSWqiC9qxKd7fI.duuYoB-JYkro9TOd_-NwXkehbCumcwyC_VNTTpNfVbw&dib_tag=se&keywords=mathematics+handbook&qid=1733585380&sprefix=mathematics+hand%2Caps%2C193&sr=8-1&linkCode=ll1&tag=engui-20&linkId=4f7b8f54bfd4a3da9ba225ca7cf7a7f8&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 712</a>).

## Continuity Equation in Cylindrical Coordinates

The continuity equation for cylindrical coordinates $(r,\theta,z)$ is

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r} \frac{\partial(\rho r u_r)}{\partial r} + \frac{1}{r} \frac{\partial(\rho u_\theta)}{\partial \theta} + \frac{\partial(\rho u_z)}{\partial z} = 0,
$$

where $u_r$, $u_\theta$, and $u_z$ are the flow velocity in the $r$-, $\theta$-, and $z$-direction respectively.

### Continuity Equation in Cylindrical Coordinates With Incompressible Flow

If the fluid is incompressible, $\rho$ is constant then the formula becomes:

$$
    \rho \left( \frac{1}{r} \frac{\partial(r u_r)}{\partial r} + \frac{1}{r} \frac{\partial(u_\theta)}{\partial \theta} + \frac{\partial u_z}{\partial z} \right) = 0
$$

 (<a href="https://www.amazon.com/Handbook-Mathematics-I-N-Bronshtein/dp/3662462206?crid=2EY5M5H7FT2IH&dib=eyJ2IjoiMSJ9.DwfzsUXXkM7e94_ArqNkD7rI9LmiHe7Nl2zScHPEeCKOU1geuDYJusAIGYH94tYi0OQ9FNiQ4yi64-rEENwjhruC-ruHygf1MlPUoylnNrhbsTbUckL03fNH60wBfKaC-Uwew6hftZIJpPJZUSt1u1-MP39iufzwGsdI5WdbGMJm9l7J_Ea3MFirPBITTnUrW7Zx80gGJl0GXL1JdcZzsvBPWQ944WSWqiC9qxKd7fI.duuYoB-JYkro9TOd_-NwXkehbCumcwyC_VNTTpNfVbw&dib_tag=se&keywords=mathematics+handbook&qid=1733585380&sprefix=mathematics+hand%2Caps%2C193&sr=8-1&linkCode=ll1&tag=engui-20&linkId=4f7b8f54bfd4a3da9ba225ca7cf7a7f8&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 712</a>).

## Continuity Equation in Spherical Coordinates

The continuity equation for spherical coordinates $(r,\theta,\varphi)$ is

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial(\rho r^2 u_r)}{\partial r} + \frac{1}{r \sin{\theta}} \frac{\partial(\rho u_\theta \sin{\theta})}{\partial \theta} + \frac{1}{r \sin{\theta}} \frac{\partial(\rho u_\varphi)}{\partial \varphi} = 0,
$$

where $u_r$, $u_\theta$, and $u_\varphi$ are the flow velocity in the $r$-, $\theta$-, and $\varphi$-direction respectively.

### Continuity Equation in Spherical Coordinates With Incompressible Flow

If the fluid is incompressible, $\rho$ is constant then the formula becomes:

$$
    \rho \left( \frac{1}{r^2} \frac{\partial (r^2 u_r)}{\partial r} + \frac{1}{r \sin{\theta}} \frac{\partial(u_\theta \sin{\theta})}{\partial \theta} + \frac{1}{r \sin{\theta}} \frac{\partial(u_\varphi)}{\partial \varphi} \right) = 0
$$

(<a href="https://www.amazon.com/Handbook-Mathematics-I-N-Bronshtein/dp/3662462206?crid=2EY5M5H7FT2IH&dib=eyJ2IjoiMSJ9.DwfzsUXXkM7e94_ArqNkD7rI9LmiHe7Nl2zScHPEeCKOU1geuDYJusAIGYH94tYi0OQ9FNiQ4yi64-rEENwjhruC-ruHygf1MlPUoylnNrhbsTbUckL03fNH60wBfKaC-Uwew6hftZIJpPJZUSt1u1-MP39iufzwGsdI5WdbGMJm9l7J_Ea3MFirPBITTnUrW7Zx80gGJl0GXL1JdcZzsvBPWQ944WSWqiC9qxKd7fI.duuYoB-JYkro9TOd_-NwXkehbCumcwyC_VNTTpNfVbw&dib_tag=se&keywords=mathematics+handbook&qid=1733585380&sprefix=mathematics+hand%2Caps%2C193&sr=8-1&linkCode=ll1&tag=engui-20&linkId=4f7b8f54bfd4a3da9ba225ca7cf7a7f8&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 712</a>).

Continuity Equation for Fluids in Different Coordinate Systems

Free vibration occurs when a spring-mounted system oscillates without external forces. Let us solve the formulas and equations of free vibration.

**Free vibration occurs when a spring-mounted system oscillates without external forces. Let us solve the formulas and equations of free vibration.**

## Undamped Free Vibration

Free vibration is the term for when a spring-mounted body is disturbed from its position of equilibrium. Undamped free vibration is ideal because all free vibration of particles, to a varying degree, dampens due to mechanical or fluid friction.

## Undamped Free Vibration Formula

Looking at the figure above and doing a free-body diagram we can deduct using Newton's second law that

$$
    m \ddot{x} + kx = 0,
$$

where $m$ is the body's mass, $x$ is the distance from equilibrium (when the force from the spring is 0), $\ddot{x}$ is the second time-derivative of $x$ and $k$ is the spring constant.

## Undamped Free Vibration Equation Solved

Looking at the ODE in the equation above we can expect an oscillatory motion so we can write $x$ as

$$
    x = A \cos{(\omega_n t)} + B \sin{(\omega_n t)}.
$$

Inserting this in the ODE we get:

$$
    m (-A \omega^2 \cos{(\omega_n t)} - B \omega^2 \sin{(\omega_n t)}) 
$$
$$
+ k (A \cos{(\omega_n t)} + B \sin{(\omega_n t)}) = 0.
$$

By solving for the sine and cosine terms we get:

$$
    -m B \omega_n^2 + k B = 0 \text{ and } -m A \omega_n^2 + k A = 0.
$$

### Natural Circular Frequency Formula

So we only have a solution if

$$
    \omega_n = \sqrt{\frac{k}{m}},
$$

which is called the natural circular frequency. Knowing that the expression for $x$ is the solution we can calculate $A$ and $B$. If we evaluate $x$ and $\dot{x}$ at $t=0$ we get $x_0 = A$ and $\dot{x} = B \omega$. Substituting these values into the solution we get

$$
    x = x_0 \cos{(\omega_n t)} + \frac{\dot{x}_0}{\omega} \sin{(\omega_n t)}.
$$

## Equilibrium Position For Vertical Spring

![Equilibrium Position For Vertical Spring](https://i.imgur.com/eP6UUaZ.png "Equilibrium Position For Vertical Spring")

If we were to flip the system mentioned above we would also have the gravitational force so now the equilibrium position will contain a nonzero spring deflection $\delta_\text{eq}$. So now the equation for Newton's second law is

$$
    m \ddot{x} + k(\delta_\text{eq} + x) - mg = 0.
$$

For the equilibrium position to be equilibrium the sum of all forces must be zero so

$$
    k \delta_\text{eq} = mg.
$$

So the equation for Newton's second law becomes the formula that we had for the horizontal case:

$$
    m \ddot{x} + kx = 0.
$$

So by defining the displacement variable to originate from the equilibrium, we can ignore the equal and opposite forces associated with equilibrium.

## Damped Free Vibration

![Damped Free Vibration](https://i.imgur.com/wSIDGZ0.png "Damped Free Vibration")

Let us introduce damping into the system in the form of a viscous damper that behaves linearly. The viscous damper consists of a cylinder filled with fluid and a piston with holes so the fluid can flow between the sides of the piston. The damping force's magnitude is proportional to the velocity of the body. The direction of the force is in the opposite direction of the velocity so we can describe the damping force on the body as $-c \dot{x}$. $c$ is called the viscous damping coefficient.

## Damped Free Vibration Formula

Doing a free-body diagram on the body gives us this equation derived from Newton's second law:

$$
    m \ddot{x} + c \dot{x} + kx = 0.
$$

### Viscous Damping Factor / Damping Ratio

We discovered that for undamped free vibration $\omega_n = \sqrt{\frac{k}{m}}$. Let us make that substitution to the equation above and let us also make this substitution:

$$
    \zeta = \frac{c}{2 m \omega_n},
$$

which is called the viscous damping factor or damping ratio which is a measurement of the severity of the damping. So now the equation becomes:

$$
    \ddot{x} + 2 \zeta \omega_n \dot{x} + \omega_n^2 x = 0.
$$

## Damped Free Vibration Equation Solved

To solve the equation above we assume a solution in the form of

$$
    x = A e^{\gamma t}
$$

and substituting this into the equation gets us

$$
    \gamma^2 + 2 \zeta \omega_n \gamma + \omega^2 = 0.
$$

This is called the characteristic equation and the roots for it is

$$
    \gamma_1 = \omega_n (-\zeta + \sqrt{\zeta^2 - 1})
$$
$$
    \gamma_2 = \omega_n (-\zeta - \sqrt{\zeta^2 - 1})
$$

Using the property of superposition due to the system's linearity we get the general solution which is the sum of the two individual solutions gained from the roots of the characteristic equation:

$$
    x = A_1 e^{\gamma_1 t} + A_2 e^{\gamma_2 t}
$$
$$
    = A_1 e^{\omega_n (-\zeta + \sqrt{\zeta^2 - 1}) t} + A_2 e^{\omega_n (-\zeta - \sqrt{\zeta^2 - 1}) t}.
$$

## Categories of Damped Motion

The damping factor $\zeta$ can have values between 0 and infinity so the expression $(\zeta^2 - 1)$ can be positive, negative, or zero. This leads to three categories of damped motion.

### Overdamped Motion

Overdamped motion is when $\zeta > 1$. This leads to the roots $\gamma_1$ and $\gamma_2$ being distinct from each other and them being real and negative. As a consequence, the equation of the solution above will approach zero when $t$ goes to infinity while not having any oscillations.

### Critically Damped Motion

When $\zeta = 1$ the motion is critically damped and the roots are equal to each other:

$$
    \gamma_1 = \gamma_2 = - \omega_n.
$$

The roots are also real and negative since $\omega$ is always positive. It can be shown that for critically damped motion the solution is given by

$$
    x = (A_1 + A_2 t) e^{-\omega_n t}.
$$

A critically damped system will also decay with no oscillations to zero when $t$ goes to infinity and it will do so faster than an overdamped system. The figure below shows the difference between the two cases. They have the same initial conditions and $\zeta = 3$ for the overdamped system.

![Overdamped vs Critically Damped](https://i.imgur.com/yRVyx2t.png "Overdamped vs Critically Damped")

### Underdamped Motion

The system is underdamped when $\zeta < 1$ so now $\zeta^2-1$ is negative and the solution becomes

$$
    x = A_1 e^{\omega_n (-\zeta + i \sqrt{1 - \zeta^2}) t} + A_2 e^{\omega_n (-\zeta - i \sqrt{1 - \zeta^2}) t}
$$
$$
    = (A_1 e^{i \omega_n \sqrt{1 - \zeta^2} t} + A_2 e^{i \omega_n \sqrt{1 - \zeta^2} t}) e^{-\zeta \omega_n t},
$$

since $\sqrt{-1} = i$. Let $\omega_d = \omega_n \sqrt{1 - \zeta^2}$ so

$$
    x = (A_1 e^{i \omega_d t} + A_2 e^{-i \omega_d t}) e^{-\zeta \omega_n t}.
$$

By using the formula $e^{\pm i x} = \cos{x} \pm i \sin{x}$ we can rewrite the previous equation as

$$
    [(A_1+A_2) \cos{(\omega_d t)} + i (A_1 - A_2) \sin{(\omega_d t)}] e^{-\zeta \omega_n t}.
$$

We define two new constants as $A_3 = (A_1 + A_2)$ and $A_3 = i (A_1 - A_2)$ so

$$
    x = [A_3 \cos{(\omega_d t)} + A_4 \sin{(\omega_d t)}] e^{-\zeta \omega_n t}.
$$

Since $\cos{(\omega_d t)}$ and $\sin{(\omega_d t)}$ has the same frequency the expression above can be written as

$$
    x = Ce^{\zeta \omega_n t} \sin{(\omega_d t + \varphi)}.
$$

To determine the constants $C$ and $\varphi$ you plug in the initial conditions for the system (for example $x(t=0) = x_0$ and $\dot{x}(t=0) = \dot{x}_0$). 

### Damped Natural Frequency and Damped Period

Analyzing our final formula above we can see it has an angular frequency of

$$
    \omega_d = \omega_n \sqrt{1 - \zeta^2}.
$$

This angular frequency is called the damped natural frequency. Using the damped natural frequency we get the damped period:

$$
    \tau_d = \frac{2 \pi}{\omega_d}.
$$

You can see an example in the figure below of an underdamped motion with $\zeta = 0.1$. The factor $C e^{-\zeta \omega_n t}$ dictates the amplitude of the motion.

![Underdamped Motion](https://i.imgur.com/ONIVGm4.png "Underdamped Motion")

### Determining The Damping Ratio With Logarithmic Decrement

We can calculate the damping ratio if we were to measure two successive amplitudes (at for example $x_1$ and $x_2$) a damped period apart:

$$
    \frac{x_1}{x_2} = \frac{C e^{-\zeta \omega_n t_1}}{C e^{-\zeta \omega_n (t_1 + \tau_d)}} = e^{\zeta \omega_n \tau_d}.
$$

Taking the natural logarithm of the previous equation we get the logarithmic decrement $\delta$:

$$
    \delta = \ln{\frac{x_1}{x_2}} = \zeta \omega_n \tau_d = \zeta \omega_n \frac{2 \pi}{\omega_n \sqrt{1 - \zeta^2}} = \frac{2 \pi \zeta}{\sqrt{1 - \zeta^2}}.
$$

Now we can solve for the damping ratio:

$$
    \zeta = \frac{\delta}{\sqrt{4 \pi^2 + \delta^2}}.
$$

Undamped and Damped Oscillations in Free Vibration

Two moving objects' motion relative to each other dictates the dynamics of their relative position, velocity, and acceleration.

**Two moving objects' motion relative to each other dictates the dynamics of their relative position, velocity, and acceleration.**

## Relative Position Formula

The equation for the relative position between an object $A$ and an object $B$ looks like this:

$$
    \vec{r}_A = \vec{r}_B + \vec{r}_{A/B},
$$

where $\vec{r}_A$ and $\vec{r}_B$ are the absolute positions of the two objects in a fixed coordinate system ($XY$ coordinate system in the figure above). $\vec{r}_{A/B}$ is the position vector of object $A$ relative to object $B$, "$A/B$" means "$A$ relative to $B$". Another way of thinking about it is that $\vec{r}_{A/B}$ is object $A$'s coordinates in a coordinate system whose origin follows object $B$ ($xy$ coordinate system in the figure above).

## Relative Velocity Formula

If we differentiate the equation for the relative position we get the velocities:

$$
    \frac{d\vec{r}_A}{dt} = \frac{d\vec{r}_B}{dt} + \frac{d\vec{r}_{A/B}}{dt},
$$

which can be rewritten as

$$
    \vec{v}_A = \vec{v}_B + \vec{v}_{A/B},
$$

where $\vec{v}_A$ and $\vec{v}_B$ are the absolute velocities of the two objects in a fixed coordinate system and $\vec{v}_{A/B}$ is the velocity of object $A$ relative to object $B$.

## Relative Velocity Due to Rotation

![Relative Velocity Due to Rotation](https://i.imgur.com/rianJ4q.png "Relative Velocity Due to Rotation")

Relative velocity can also be used between two points on a rotating rigid body, see the figure above. This means that the distance ($\vec{r}_{A/B}$) between the two points remains constant. Thus the relative velocity vector $\vec{v}_{A/B}$ depends on the angular velocity $\omega$ around the origin of the translating coordinate system attached to point $B$ according to the following formula:

$$
    v_{A/B} = \omega r_{A/B}.
$$

This can be written in vector form with this:

$$
    \vec{v}_{A/B} = \vec{\omega} \times \vec{r}_{A/B},
$$

where $\vec{\omega}$ is the angular velocity vector normal to the plane of motion. Therefore the relative velocity is always perpendicular to $\vec{r}_{A/B}$. So the formula for the velocity of point $A$ is:

$$
    \vec{v}_A = \vec{v}_B + \vec{\omega} \times \vec{r}_{A/B}.
$$

## Relative Acceleration Formula

If we differentiate the equation for relative position a second time we get the accelerations:

$$
    \frac{d^2\vec{r}_A}{dt^2} = \frac{d^2\vec{r}_B}{dt^2} + \frac{d^2\vec{r}_{A/B}}{dt^2},
$$

which can be rewritten as

$$
    \vec{a}_A = \vec{a}_B + \vec{a}_{A/B},
$$

where $\vec{a}_A$ and $\vec{a}_B$ are the absolute accelerations of the two objects in a fixed coordinate system and $\vec{a}_{A/B}$ is the acceleration of object $A$ relative to object $B$. We can note here that if object $B$ is moving at a constant velocity, the acceleration of object $A$ will be the same in both the fixed coordinate system and the translating coordinate system attached to object $B$.

## Relative Acceleration Due to Rotation

![Relative Acceleration Due to Rotation](https://i.imgur.com/KbdBs8w.png "Relative Acceleration Due to Rotation")

As mentioned for the relative velocity due to rotation the distance ($\vec{r}_{A/B}$) between the two points remains constant, as can be seen in the figure above. This means that the observer moving with point $B$ perceives point $A$ to have circular motion around point $B$. This means the relative acceleration will have a normal component parallel with $\vec{r}_{A/B}$ from $A$ pointing to $B$. This is because of the change in direction of $\vec{r}_{A/B}$, like centripetal acceleration, so it has the following expression:

$$
    (a_{A/B})_n = \frac{v_{A/B}^2}{r_{A/B}} = \omega^2 r_{A/B}
$$

It will also have a tangential component due to the change in magnitude of $\vec{r}_{A/B}$ and it looks like this:

$$
    (a_{A/B})_t = \frac{dv_{A/B}}{dt} = \alpha r_{A/B},
$$

where $\alpha$ is the angular acceleration. In vector form both equations look like this:

$$
    (\vec{a}_{A/B})_n = \vec{\omega} \times (\vec{\omega} \times \vec{r}_{A/B})
$$
$$
    (\vec{a}_{A/B})_t = \vec{\alpha} \times \vec{r}_{A/B}.
$$

So the formula for the relative acceleration due to rotation looks like this:

$$
    \vec{a}_A = \vec{a}_B + \vec{\omega} \times (\vec{\omega} \times \vec{r}_{A/B}) + \vec{\alpha} \times \vec{r}_{A/B}.
$$

Relative Motion, Velocity, and Acceleration

The selection rules are important for classifying beta decay's transitions. They are based on the conservation of angular momentum and parity during decay.

**The selection rules are important for classifying beta decay's transitions. They are based on the conservation of angular momentum and parity during decay.**

## Beta + vs Beta - Decay

There exist different kinds of beta decay, two of them being beta + ($\beta^+$) and beta - ($\beta^-$) decay. In $\beta^+$ decay, a proton $p^+$ decays into a neutron $n$, a positron $e^+$, and an electron neutrino ${\nu}_e$:

$$
    p^+ \rightarrow n + e^+ + \nu_e.
$$

In $\beta^-$ decay, a neutron $n$ decays into a proton $p^+$, an electron $e^-$, and an electron antineutrino $\bar{\nu}_e$:

$$
    n \rightarrow p + e^- + \bar{\nu}_e
$$

(<a href="https://www.amazon.com/Introductory-Nuclear-Physics-Kenneth-Krane/dp/047180553X?crid=3DKB5L641CP1H&dib=eyJ2IjoiMSJ9.b55RMLFLu7cDCmPl8nhU8y5W_-6HlpX2wDw39nPZ1gaN4RAfZ76q9lMnPBmCIjo4zQXEAwv2euBUQbpnJZVeJPfDN0LeUrVgxNgXMty__xENinxsg8qo5rJKivcKmdI1lu-i4El8qaDe9zcI59eKTw.6B2FmeVbIp-7y0QDWMhZl0n8IsuynNyOvQA9spg-pjw&dib_tag=se&keywords=introductory+nuclear+physics+krane&qid=1733586348&sprefix=introductory+nuclear+p%2Caps%2C220&sr=8-1&linkCode=ll1&tag=engui-20&linkId=7439c3f2582231001d017c8a518b505d&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 272</a>).

## Fermi and Gamow-Teller Transitions

For $\beta$ decay, the total spin quantum number $S$ is obtained by summing the spins of the emitted particles. Since both the electron and the neutrino have a spin of 1/2, the sum will therefore be either 0 (antiparallel) or 1 (parallel), depending on the orientation of the spins. $S=0$ is called a Fermi transition, while $S=1$ is called a Gamow-Teller transition. 

If we take a look at the change in angular momentum of the nucleus

$$
    \Delta I = |I_i - I_f|,
$$

where $I_i$ is the angular momentum of the initial nuclear state and $I_f$ is the angular momentum of the final nuclear state. We see that $\Delta I$ does not change ($\Delta I = 0$) for a Fermi transition since $S=0$. When the electron and neutrino spins are parallel (Gamow-Teller decay) they carry a total angular momentum of 1 unit so the conversation of angular looks like this in vector form:

$$
    \vec{I_i} = \vec{I_f} + \vec{1}.
$$

This is only possible if $\Delta I = 0$ or $1$. However, it is not possible when $I_i = I_f = 0$ ($0 \rightarrow 0$), in that case only the Fermi transition can contribute (<a href="https://www.amazon.com/Introductory-Nuclear-Physics-Kenneth-Krane/dp/047180553X?crid=3DKB5L641CP1H&dib=eyJ2IjoiMSJ9.b55RMLFLu7cDCmPl8nhU8y5W_-6HlpX2wDw39nPZ1gaN4RAfZ76q9lMnPBmCIjo4zQXEAwv2euBUQbpnJZVeJPfDN0LeUrVgxNgXMty__xENinxsg8qo5rJKivcKmdI1lu-i4El8qaDe9zcI59eKTw.6B2FmeVbIp-7y0QDWMhZl0n8IsuynNyOvQA9spg-pjw&dib_tag=se&keywords=introductory+nuclear+physics+krane&qid=1733586348&sprefix=introductory+nuclear+p%2Caps%2C220&sr=8-1&linkCode=ll1&tag=engui-20&linkId=7439c3f2582231001d017c8a518b505d&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 289-292</a>).

## Allowed and Forbidden Transitions

The orbital angular momentum $L_\beta$ can take on all positive integers and also dictates the type of transition. When $L_\beta=0$ the decay is called allowed, while the remaining values are called forbidden (because they are generally less probable than allowed decays), with $L_\beta=1$ being the first-forbidden decay, $L_\beta=2$ being the second-forbidden, and so on. The allowed and all the forbidden decays have each a Fermi and Gamow-Teller type, and the allowed Fermi decay is sometimes called superallowed (<a href="https://www.amazon.com/Introductory-Nuclear-Physics-Kenneth-Krane/dp/047180553X?crid=3DKB5L641CP1H&dib=eyJ2IjoiMSJ9.b55RMLFLu7cDCmPl8nhU8y5W_-6HlpX2wDw39nPZ1gaN4RAfZ76q9lMnPBmCIjo4zQXEAwv2euBUQbpnJZVeJPfDN0LeUrVgxNgXMty__xENinxsg8qo5rJKivcKmdI1lu-i4El8qaDe9zcI59eKTw.6B2FmeVbIp-7y0QDWMhZl0n8IsuynNyOvQA9spg-pjw&dib_tag=se&keywords=introductory+nuclear+physics+krane&qid=1733586348&sprefix=introductory+nuclear+p%2Caps%2C220&sr=8-1&linkCode=ll1&tag=engui-20&linkId=7439c3f2582231001d017c8a518b505d&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 289-292</a>). 

## Parity in Beta Decay

$L_\beta$ also determines the parity change of the wave function for beta decay according to $\Delta \pi=(-1)^{L_\beta}$. So for example, allowed transitions' parity does not change.

## First-forbidden Decay

First-forbidden decays, $L_\beta = 1$, like allowed decays, have both a Fermi and Gamow-Teller type. For the Fermi transition type, we couple $S=0$ with $L_\beta = 1$ to get a total of 1 unit of angular momentum carried by the beta decay so the change in angular momentum becomes $\Delta I = 0$ or $1$ (but no $0 \rightarrow 0$ transition). For the Gamow-Teller type $S=1$ combined with $L_\beta = 1$ gives a total of 0, 1, or 2 units of total angular momentum so $\Delta I = 0$, $1$, or $2$. Since $L_\beta = 1$ is an odd number the parity changes so $\Delta \pi=(-1)^{L_\beta} = $ yes (<a href="https://www.amazon.com/Introductory-Nuclear-Physics-Kenneth-Krane/dp/047180553X?crid=3DKB5L641CP1H&dib=eyJ2IjoiMSJ9.b55RMLFLu7cDCmPl8nhU8y5W_-6HlpX2wDw39nPZ1gaN4RAfZ76q9lMnPBmCIjo4zQXEAwv2euBUQbpnJZVeJPfDN0LeUrVgxNgXMty__xENinxsg8qo5rJKivcKmdI1lu-i4El8qaDe9zcI59eKTw.6B2FmeVbIp-7y0QDWMhZl0n8IsuynNyOvQA9spg-pjw&dib_tag=se&keywords=introductory+nuclear+physics+krane&qid=1733586348&sprefix=introductory+nuclear+p%2Caps%2C220&sr=8-1&linkCode=ll1&tag=engui-20&linkId=7439c3f2582231001d017c8a518b505d&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 291</a>).

## Second-forbidden Decay

Using a similar methodology as the first-forbidden decay if we couple $L_\beta = 2$ and $S=0$ we in principle get $\Delta I = 0$, $1$, or $2$ for the Fermi type, and for the Gamow-Teller type ($L_\beta = 2$ and $S=1$) we in principle get $\Delta I = 0$, $1$, $2$, or $3$. However, since the parity does not change ($\Delta \pi=(-1)^{L_\beta} = $ no) the $\Delta I = 0$ and $1$ cases fall within the allowed decay's selection rules so we can exclude them since the allowed variations are much more probable to occur (<a href="https://www.amazon.com/Introductory-Nuclear-Physics-Kenneth-Krane/dp/047180553X?crid=3DKB5L641CP1H&dib=eyJ2IjoiMSJ9.b55RMLFLu7cDCmPl8nhU8y5W_-6HlpX2wDw39nPZ1gaN4RAfZ76q9lMnPBmCIjo4zQXEAwv2euBUQbpnJZVeJPfDN0LeUrVgxNgXMty__xENinxsg8qo5rJKivcKmdI1lu-i4El8qaDe9zcI59eKTw.6B2FmeVbIp-7y0QDWMhZl0n8IsuynNyOvQA9spg-pjw&dib_tag=se&keywords=introductory+nuclear+physics+krane&qid=1733586348&sprefix=introductory+nuclear+p%2Caps%2C220&sr=8-1&linkCode=ll1&tag=engui-20&linkId=7439c3f2582231001d017c8a518b505d&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 292</a>). 

## Beta Decay Transition Rules

A table consolidating all the transition rules for beta decay mentioned above is shown below.

![Beta Decay Transitions, Allowed, Forbidden, Fermi, Gamow-Teller, Parity Change](https://i.imgur.com/CPhtzE3.png "Beta Decay Transitions, Allowed, Forbidden, Fermi, Gamow-Teller, Parity Change")

Allowed, Forbidden, Fermi, and Gamow-Teller Beta Decay Transition

What is the smallest radius a rolling disk can rotate around its precession center? Let's explore the disk's angular velocity and precession rate.

**At how small radius can a rolling disk rotate around its precession center? What is the disk's angular velocity and its precession rate? Let us look at the physics that determines these values.**

## Given Variables

We assume to know the values of the following variables:

- $R$: The distance from the center of precession to the contact point of the disk
- $r$: The radius of the disk
- $g$: The acceleration due to gravity
- $\theta$: The angle between the disk and the vertical to the ground

## Variables, Units, and Assumptions for Searched Values

We define the following assumptions and variables which we want to calculate:

- $\omega$: The angular velocity (unit rad/s) around the center of mass of the disk
- $\Omega$: The precession rate (unit rad/s), also known as the precession angular velocity, of the disk
- The minimum value of $R$ for the phenomenon to occur
- Lastly, we assume that the disk is homogeneous.

## Relation Between Disk's Angular Velocity and Precession Rate

The speed around the center of mass of the disk, $\omega r$, is the same as that around the precession point, $\Omega R$. This leads to the following:

$$
    \omega=\frac{\Omega R}{r}.
$$

## Equation For Precession Rate

Projecting the force of gravity that acts upon the disk on the local coordinate system shown in the figure above the moment around the contact point becomes the following:

$$
    \vec{M}= r \hat{z} \times mg (-\sin{\theta} \hat{x} - \cos{\theta} \hat{z}) = -mgr \sin{\theta} \hat{y}.
$$

The center of mass is closer to the center of precession than the contact point so the the velocity of the center of mass is

$$
    \vec{v} = \Omega (R - r \sin{\theta}) \hat{y}.
$$

The rotational velocity around the center of mass in the vector in the local coordinate system is

$$
    \vec{\omega} = \begin{bmatrix}
    \frac{\Omega R}{r} \\
    0 \\
    0
    \end{bmatrix}
    .
$$

The precession rate translated to the local coordinate system:

$$
    \vec{\Omega} = \begin{bmatrix}
    \Omega \sin{\theta} \\
    0 \\
    \Omega \cos{\theta}
    \end{bmatrix}
    .
$$

We also need the [moment of inertia around the center of mass of the disk](https://www.engineeringguidance.com/the-moment-of-inertia-of-a-rod-and-a-disk), which is

$$
    I = \begin{bmatrix}
    2mr^2k & 0 & 0 \\
    0 & mr^2k & 0 \\
    0 & 0 & mr^2k
    \end{bmatrix},
$$

where $k=\frac{1}{4}$. We can now use expressions in vector form so the angular momentum around the contact point becomes

$$
    \vec{H} = I (\vec{\omega} + \vec{\Omega}) + r \hat{z} \times \vec{v} 
$$
$$
= mr \Omega ((r \sin{\theta} - R)(2k +2) \hat{x} + r(k+1) \cos{\theta} \hat{z}).
$$

Then we take the time derivative of the angular momentum at the contact point:

$$
        \frac{d \vec{H}}{dt} = \frac{\partial \vec{H}}{\partial t} + \vec{\Omega} \times \vec{H}
$$

Since the components are not time-dependent, the time-related derivative $\frac{\partial \vec{H}}{\partial t}$ is zero. Therefore the derivative of the angular momentum around the contact point becomes

$$
        \vec{\Omega} \times \vec{H} 
$$
$$
= -\Omega^2 r m \cos{\theta} (rk \sin{\theta} + (R - r \sin{\theta})(2k + 1)) \hat{y}.
$$

The derivative of the angular momentum at the contact point is equal to the moment at the contact point so

$$
        -mgr \sin{\theta} \hat{y} 
$$
$$
=  -\Omega^2 r m \cos{\theta} (rk \sin{\theta} + (R - r \sin{\theta})(2k + 1)) \hat{y}.
$$

With that equation, $\Omega$ can be solved and we get the equation for the precession rate:

$$
    \Omega = \sqrt{\frac{g \tan{\theta}}{rk \sin{\theta} + (R - r \sin{\theta})(2k + 1)}}.
$$

## Minimum value of $R$

When the equation above becomes undefined, that is, when the denominator becomes zero, the minimum value of $R$ for precession to occur. So we have have:

$$
    0 = rk \sin{\theta} + (R - r \sin{\theta})(2k + 1).
$$

So if we solve for $R$ we get:

$$
    R = \frac{5}{6} r \sin{\theta}.
$$

## Analysis of Solution

During these calculations, an assumption regarding friction has been made. Some friction must exist to keep the disk in a circular path. Furthermore, these calculations also assume that the disk has no height, which means that $I_{22}$ and $I_{33}$ in $I$ do not have a height term. On top of that, the edge radius in these calculations cannot exist for the length vectors to be correct and at the same time, a larger edge radius would provide a larger contact area, which would lead to greater maximum friction.

If $\theta = 0$, no moment would arise since the gravitational and normal forces would align and precession could not occur. If the velocity were also very low, only a slight instability would cause a moment to arise and since the time derivative of the angular momentum must equal the moment, this means that a large change occurs in the angular momentum when it is small. This would likely lead to the disk falling.

Precession Rate of a Rolling Disk

A charge density is spherically symmetric when it depends only on distance from a point, not direction. We’ll walk through finding 40Ca's nucleus charge density.

**A charge density is considered spherically symmetric when its value depends only on the distance from a specific point space and not the direction. Consequently, angles do not matter, so it looks the same if the system is rotated. In this post, we will walk through a method of determining $^{40}$Ca's nucleus charge density.**

## Charge Density of 40Ca

As an example, let us determine the spherically symmetric charge density

$$
    \rho_\text{ch}(r;\mathbf{X}) = \frac{\rho^0_{\text{ch}}}{1+e^{(r-b)/a}}
$$

for $^{40}$Ca. Therefore, the constants $\mathbf{X} = [\rho^0_{\text{ch}},a,b]$ needs to be determined. 

## Doubly Magic Nuclei

$^{40}$Ca has a double magic number, meaning the atomic nucleus is more stable than other isotopes. Magic atomic nuclei have very small electric quadrupole moments, making them more or less spherical, which is their most tightly bound state.

## Measuring 40Ca's Charge Density

This is done by elastically irradiating the target $^{40}$Ca with electrons having an energy of 250 MeV. The target is assumed to be stationary, and no recoil is assumed to occur since the target is much heavier than the electron. Since the target is assumed to be stationary in the LAB frame and no recoil is assumed for the Ca atom, the COM (center-of-momentum) frame does not move in the LAB frame. Thus, energies, scattering cross-sections, and angles are approximately the same in both frames.

Differential cross-sections $\mathcal{E}_i$ has been measured for the electrons to be scattered at an angle $\theta$. These cross-sections have measurement errors $\delta_i$. At the same time, the cross-section can be theoretically expressed as

$$
    \mathcal{T}(\mathbf{X}) = \left(\frac{\sigma}{\Omega} \right)_\text{Mott} |F|^2 
$$
$$
= \frac{Z_p^2 Z_t^2 \alpha^2 (\hbar c)^2}{4 \beta^4 E^2 \sin^4(\theta/2)} [1-\beta^2 \sin^2(\theta/2)] \cdot |F|^2,
$$

where $F$ is the form factor, $\alpha$ is the fine structure constant, $\beta = v/c$ with $v$ being the electron's speed, $E$ is the electron's energy, $Z_p$ and $Z_t$ are the charges of the electron and the target, respectively, and $\theta$ is the scattering angle. In order for the measured cross-sections to agree with the theoretical values, minimization was performed on

$$
    \mathcal{X} ^2(\mathbf{X}) = \sum_{i=1}^{N_\theta} \left(\frac{\mathcal{T}(\mathbf{X})-\mathcal{E}_i}{\delta_i} \right)^2.
$$

This needs to be done with the constraint

$$
    Z_t = \frac{4 \pi}{e} \int_0^\infty r^2 \rho_\text{ch}(r;\mathbf{X}) dr.
$$

The constraint ensures that the total charge of the charge density equals the charge of the atomic nucleus $Z_t e$. Without the constraint, the form factor and charge radius would yield incorrect values. The charge density could potentially also become negative, which is not possible. This was done using **scipy**'s **basinhopping**, where a starting point

$$
    \mathbf{X} = [0.06 \text{ e fm}^{-3}, 1.07 \cdot 40^{1/3} \text{ fm}, 0.54 \text{ fm}]
$$

was used.

## Measurements vs Theoretical Value

The vector $\mathbf{X}_*$ that gives the minimum $\mathcal{X}_* ^2$ sum is $\mathbf{X} = [0.077 \text{ e fm}^{-3}, 3.8 \text{ fm}, 0.47 \text{ fm}]$. In the figure below, the charge density can be observed. In the figure below that, a plot of the theoretical and measured cross-sections can also be seen.

![Approximate charge density](https://i.imgur.com/Zr5Vqo0.png "Approximate charge density")

![Cross-section measured vs theoretical](https://i.imgur.com/handYRg.png "Cross-section measured vs theoretical")

## 40Ca's rms charge radius

Using the charge function, the rms charge radius can be calculated as follows:

$$
    \sqrt{\langle r^2 \rangle} = \sqrt{\frac{4 \pi}{Z_t e} \int_0^\infty r^4 \rho_\text{ch}(r;\mathbf{X}_*) dr} \approx 3.41 \text{ fm}.
$$

This differs from the accepted value of 3.4776(19) fm. The difference may be due to the fact that the charge density function is an approximation. The "true" charge distribution has a similar form, but it does not precisely follow the model used. Furthermore, there is uncertainty in the value of $\rho^0_{\text{ch}}$ at small $r$.

Spherical Symmetric Charge Distribution of Nucleis

The ARMA process, a time series model, is key in forecasting. We'll explore its definition, stationarity, causality, invertibility, and model order.

**The ARMA process, short for autoregressive moving average, is one the most basic time series models but it is an important building block in forecasting. Let us look at its definition and some of its properties such as stationarity, causality, and invertibility but also how to choose which order your ARMA model should have.**

## ARMA Process Definition

The ARMA($p,q$) process $X_t$ has the following defintion:

$$
    X_t - \sum_{i=1}^q \phi_j X_{t-i} = Z_t + \sum_{i=1}^q \theta_j Z_{t-i},
$$

where $Z_t \sim \text{WN}(0,\sigma^2)$ (white noise with mean 0 and variance $\sigma^2$) and $1 -\sum_{i=1}^q \phi_i z^i$ and $1 + \sum_{i=1}^q \theta_i z^i$ have no common zeros. You might also define the ARMA model in the following short form:

$$
    \phi(B) X_t = \theta(B) Z_t,
$$

where $B$ is the [lag operator](https://www.engineeringguidance.com/time-series-seasonality-and-trend) defined as $B X_t = X_{t-1}$. The two functions are defined as:

$$
    \phi(z) = \sum_{i=1}^q \phi_i z^i
$$

and

$$
    \theta (z) = \sum_{i=1}^q \theta_i z^i
$$

## When is an ARMA Process Stationary?

An ARMA process $X_t$ has a unique [(weakly) stationary](https://www.engineeringguidance.com/time-series-white-noise) if and only if $\phi(z) \neq 0$ for all $z \in \mathbb{C}$ when $|z| = 1$

## When is an ARMA Process Causal?

An ARMA process $X_t$ is causal if there exist a real value sequence $(\varepsilon_i, i \in \mathbb{N}_0)$ such that for all $t \in \mathbb{Z}$:

$$
    \sum_{i=0}^\infty |\varepsilon_i| < \infty
$$

and

$$
    X_t = \sum_{i=0}^\infty \varepsilon_i Z_{t-i}.
$$

In other words, $X_t$ can be described as an MA$(\infty)$ process. Equivalently an ARMA process $X_t$ is causal if and only if $\phi(z) \neq 0$ for all $z \in \mathbb{C}$ when $|z| \leq 1$.

## When is an ARMA Process Invertible?

An ARMA process $X_t$ is invertible if there exist a real value sequence $(\pi_i, i \in \mathbb{N}_0)$ such that for all $t \in \mathbb{Z}$:

$$
    \sum_{i=0}^\infty |\pi_i| < \infty
$$

and

$$
    Z_t = \sum_{i=0}^\infty \pi_i X_{t-i}.
$$

In other words, $X_t$ can be described as an AR$(\infty)$ process. Equivalently an ARMA process $X_t$ is invertible if and only if $\theta(z) \neq 0$ for all $z \in \mathbb{C}$ when $|z| \leq 1$.

## Choosing The Order of The ARMA Process

It’s always possible to fit an ARMA$(p,q)$ model with excessively large $p$ and $q$ values, but this isn’t advantageous for forecasting. While this approach often yields a small estimated white noise variance, the mean squared error of forecasts is also affected by errors in parameter estimation. To address this, we introduce a "penalty factor" to discourage the selection of overly complex models.

## AICC Criterion

We start by introducing the AICC criterion, where AIC stands for Akaike’s Information Criterion and the last C denotes biased-corrected. The AICC estimates the Kullback-Leibler divergence, which measures how different one probability distribution is from another, comparing the estimated distribution to the actual data distribution. It assumes that $Z \sim \text{IID} \mathcal{N} (0, \sigma^2)$, but it shows robustness to moderate deviations from normality, such as when $Z_t$ follows a $t$-distribution.

The AICC criterion says that you have to choose $p$, $q$, $\phi_p$ and $\theta_q$ to minimize the following:

$$
    -2 \ln{L(\phi_p, \theta_q, S(\phi_p,\theta_q)/n) + 2n \frac{p+q+1}{n-p-q-2}},
$$

where $\phi_p = (\phi_1, ... , \phi_p)$ and $\theta_p = (\theta_1, ... , \theta_p)$. To apply this criterion in practice, we fit a wide range of models with varying orders (p, q) to the data and select the one that minimizes the negative log-likelihood, adjusted by the penalty factor $2n \frac{p+q+1}{n-p-q-2}$

## BIC Criterion

One issue with the AICC criterion is that the estimators for p and q aren’t consistent; they don’t converge almost surely to the true values. In contrast, consistent estimators can be derived using the Bayesian Information Criterion (BIC), which also penalizes the selection of large p and q values, helping to prevent overfitting.

The BIC criterion says that you have to choose $p$ and $q$ to minimize the following:

$$
    (n-p-q) \ln{\frac{n\hat{\sigma}^2}{n-p-q}}+n (1+\ln{\sqrt{2 \pi}}) 
$$
$$
+ (p+q) \ln{\frac{\sum_{t=1}^n X_t^2 - n \hat{\sigma}^2}{p-q}},
$$

where $\hat{\sigma}^2$ is the maximum likelihood estimate of the white noise variance. A downside of the BIC is its efficiency in finding minimizers. While selecting the model that minimizes the AICC is asymptotically efficient for causal and invertible ARMA processes, this isn’t the case for the BIC. In this context, efficiency means that minimizing the AICC will, in the long run, lead to a model with the lowest one-step-ahead prediction errors.

When using constrained maximum likelihood estimation—where certain coefficients are assumed to be zero during the estimation process—the term $p+q+1$ is replaced by $m$, which represents the number of non-zero coefficients, in both the AICC and BIC.

The ARMA Process - Stationary, Causal, and Invertible

S-parameters measure a network's response at its ports, capturing reflected and transmitted waves, helping determine its behavior in new environments.

**S-parameters, short for scattering parameters, are a useful and easy way to measure a network's response to input at its terminals (ports). With the S-parameters of a network determined, the behavior of the network in a new external environment can be determined. What is measured is the reflected and transmitted waves going through the network's ports.**

## S-parameter Formula

To calculate the S-parameters of a network one can use the normalized incident complex voltage wave $a_i$ at the $i$:th port defined as

$$
    a_i = \frac{V_i + Z_i I_i}{2 \sqrt{|\text{Re} Z_i|}}
$$

and the normalized reflected complex voltage wave $b_i$ defined as

$$
    b_i = \frac{V_i - \overline{Z_i} I_i}{2 \sqrt{|\text{Re} Z_i|}},
$$

where $V_i$ is the terminal voltage, $I_i$ the terminal current and $Z_i$ is an arbitrary reference impedance (and $\overline{Z_i}$ is it's complex conjugate). Then the S-parameter $S_{ij}$ is defined as

$$
    S_{ij} = \frac{b_i}{a_j}.
$$

## S-parameters For a 2-port Network

For a 2-port network, seen in the figure above, there are four different S-parameters. These can be used to describe the relation between incident waves $a_1$ and $a_2$ and the reflected waves $b_1$ and $b_2$ in the following

$$
    \left\{
    \begin{array}{cc}
         b_1 = & S_{11} a_1 + S_{12} a_2 \\
         b_2 = & S_{21} a_1 + S_{22} a_2
    \end{array}
    \right.
    .
$$

To get $S_{11}$ the second port has to be terminated with a matched load $Z_L = Z_0$ and $a_2 = 0$ so

$$
    S_{11} = \left. \frac{b_1}{a_1} \right|_{a_2 = 0}.
$$

To get $S_{22}$ the first port has to be terminated with a matched load $Z_S = Z_0$ and $V_S = 0$ so

$$
    S_{22} = \left. \frac{b_2}{a_2} \right|_{a_1 = 0}.
$$

To get $S_{21}$ the second port has to be terminated with a matched load $Z_L = Z_0$ and $a_2 = 0$ so

$$
    S_{21} = \left. \frac{b_2}{a_1} \right|_{a_2 = 0}.
$$

To get $S_{12}$ the first port has to be terminated with a matched load $Z_S = Z_0$ and $a_1 = 0$ so

$$
    S_{12} = \left. \frac{b_1}{a_2} \right|_{a_1 = 0}.
$$

Using these formulas and by choosing a real reference impedance $Z_0$ (so $Z_0=\overline{Z_0}$ we get

$$
    S_{11} = \frac{b_1}{a_1} = \frac{\frac{V_1}{I_1} - Z_0}{\frac{V_1}{I_1} + Z_0} = \frac{Z_1 - Z_0}{Z_1 + Z_0}.
$$

Therefore we also get

$$
    Z_1 = Z_0 \frac{1 + S_{11}}{1 - S_{11}},
$$

where at port 1, $Z_1 = \frac{V_1}{I_1}$ is the input impedance.

## S-parameters For Antennas

![S-parameters For Antennas](https://i.imgur.com/sRWuN8i.png "S-parameters For Antennas")

S-parameters for an $N$:th element antenna array is calculated by

$$
    S_{mn} = \left. \frac{V_m^-}{V_n^+} \right|_{V_k^+ = 0 \text{ for } k \neq n}
$$

where $V_m^+$ are the incident and $V_m^-$ the reflected voltage wave amplitudes for the $m$:th array element seen in the figure above.

S-parameters in Circuits and Antennas

In this post, vehicle dynamics parameters, such as slip, pitch, roll, and ride height are described which can affect the aerodynamic performance of a car.

**In this post, vehicle dynamics parameters, such as slip, pitch, roll, and ride height are described which can affect the aerodynamic performance of a car.**

## Steering Angle

The steering angle is the angle between the direction of a wheel and the car's longitudinal axis. The steering angle for each front wheel can be different if the car has Ackermann steering or anti-Ackermann steering. An illustration can be found in the figure above.

## Body Slip Angle

![Body slip angle](https://i.imgur.com/rZFkoDP.png "Body slip angle")

The body slip angle is a yaw angle meaning it is the angle that is created when the car rotates around its vertical axis. The rotational axis goes through the car's center of gravity. The body slip angle is the angle between the velocity vector of the car and the car's longitudinal direction. A visual representation of the body slip angle can be seen in the figure above.

## Pitch Angle

![Pitch angle](https://i.imgur.com/gckbgDJ.png "Pitch angle")

The pitch angle is the angle created if the vehicle's sprung mass (everything but the wheels) rotates around its lateral axis, as shown in the figure above. The vehicle rotates around what is called its pitch center.

## Roll Angle

![Roll angle](https://i.imgur.com/RsVYKF6.png "Roll angle")

The roll angle is the angle that is created if the vehicle's sprung mass rotates around its longitudinal axis, see the figure above. The vehicle rotates around what is called its roll center.

## Ride Height

![Ride height](https://i.imgur.com/iGjAZla.png "Ride height")

The ride height is the smallest distance between the ground and the car's floor when it is stationary, see the figure above.
    
## Aerodynamic Influence on Vehicle Dynamics

Vehicle dynamics is a large and complex subject. Mass distribution, tires, suspension setup, and aerodynamics are examples of areas that greatly affect the dynamics of a car.

One of the main aims of aerodynamics on a race car is to increase its negative lift force (downforce). Greater downforce makes it possible for the tires to produce larger forces (grip) while accelerating/braking (longitudinal forces) and cornering (lateral forces). The basic formula describing this relationship is $F_{\text{long/lat}} = F_z \mu$, where $F_{\text{long/lat}}$ is the longitudinal or latitudinal force, $F_z$ the normal force which is increased with increased downforce and $\mu$ the friction coefficient. Consequently, the speed of the car can be increased in most driving scenarios or the same grip can be achieved with less tire slip, i.e. smaller slip angle, when compared to a car with less downforce. This helps reduce heat production and therefore preserve tire life. Furthermore, smaller slip angles are easier to control for the driver. The car's stability in this case is improved compared to larger slip angles and aerodynamic parts tend to work more efficiently, improving the balance.

In most cases, a weight and aerodynamic distribution close to 50/50 on the car's rear and front axle is desired. A car with the center of gravity (COG) closer to the front axle usually generates more understeer. Conversely, the aerodynamic center of pressure (COP or aero balance), the point at which the resultant downforce is located, closer to the front usually generates more oversteer. A car with a front-biased COG can therefore be counterbalanced by placing the COP forward. However, it is favorable in some cases to design the aerodynamic package with a COP slightly rear-biased. This can generate understeer, but understeer is considered more stable and easier to control when compared to oversteer.

## The Front Wing

Vehicle dynamics parameters described above are closely related to the aerodynamic package on a given car. While the generated downforce affects body slip, pitch, and roll angles, these parameters in turn affect the aerodynamic efficiency. Pitch angles for example can have a large impact on vehicle balance. When braking and pitching forwards, the front wing's/front splitters's downforce is generally increased. 

However, other factors such as already low static ride height or an extreme pitch angle might cause a decrease in downforce for the front wing. This also results in less airflow to the car's floor and diffuser, decreasing the middle and rear downforce. If the front wing downforce increases, the aerodynamic balance shifts forward, decreasing the stability of the car. While this is considered unfavorable, increased grip on the front axle also results in improved braking capabilities. On the contrary, if the downforce on the front wing decreases the balance shifts rearwards and the opposite occurs. It is more difficult to describe the influence of steering and roll angles. These affect the aerodynamics in an asymmetric manner and can cause unpredictable behavior in the car. The same is true for body slip angles, or yaw angles, where the airflow hits the car with an angle from the longitudinal axis.

## The Diffuser

A diffuser expands the airflow underneath the car's floor. This expansion forces the flow upstream to increase in speed in order to fill up the larger volume in the diffuser, decreasing the pressure upstream. In motorsport, it has been observed that diffusers are particularly effective aero devices. The reason is that diffusers increase the vehicle's downforce under the whole floor. Downforce in this region is close to the center of the car and distributed over a large area. This results in a good balance contribution and makes the vehicle less sensitive to pitch or other phenomena that affect the balance. 

In order for a diffuser to work efficiently, the front wing needs to supply the floor and diffuser area with sufficient airflow. Optimally the airflow's velocity is large and laminar. Another option for increasing the airflow is by feeding the diffuser with air from the sides of the vehicle. This can however interfere with the low-pressure diffuser region and increase the pressure, resulting in a decrease of downforce. A compromise between keeping the floor sealed and increasing the airflow is usually optimal. This airflow can also be used to induce vortices.

Relation Between Vehicle Dynamics and Aerodynamics

Signal distortion is unwanted changes in the output, often harmonic, from slew rate or clipping, occurring at multiples of the fundamental frequency.

**Signal distortion is unwanted changes in the output signal relative to the input signal. Some of these distortions can be harmonic, for example, from slew rate or clipping, meaning that disturbances occur at multiples of the fundamental frequency.**

## Total Harmonic Distortion

Distortion can be measured in various ways, but a common measurement is total harmonic distortion (THD). THD is defined as

$$
    \text{THD} = \frac{\sqrt{V_2^2+V_3^2+ ... +V_n^2}}{V_1},
$$

where $V_n$ is the n-th harmonic amplitude in the frequency domain. $V_1$ is called the fundamental frequency. Generally, only the first five or six harmonic amplitudes are significant for THD, and often the first three are sufficient in practice without greatly impacting the margin of error.

## Slew Rate of Operational Amplifier

Harmonic distortion manifests in various forms. One example is the slew rate of an operational amplifier. The slew rate denotes the maximum rate of change of the output signal, dictated by the specified slew rate value of the operational amplifier. When this rate of change is exceeded, the waveform is replaced by a linear ramp that adheres to the maximum slew rate until it eventually intersects with the original waveform again, as illustrated in the figure above. This phenomenon can lead to significant alterations in signal fidelity, particularly in high-frequency applications, where the ability to accurately track rapid changes in the input signal is critical for maintaining overall performance.

## Clipping in Operational Amplifier

Another form of harmonic distortion is clipping. Clipping can be categorized into soft clipping and hard clipping. 

### Hard Clipping

![Hard Clipping](https://i.imgur.com/GlgkIKr.png "Hard Clipping")

Hard clipping means that the output signal can only reach a certain amplitude, resulting in the waveform appearing "chopped off." For a negatively feedbacked operational amplifier, this amplitude depends, among other factors, on the power supply voltage. The output voltage cannot exceed the operational amplifier's supply voltage. Another way to achieve hard clipping is by using diodes between the output voltage and ground, as the output voltage in that case cannot exceed the threshold voltage

### Soft Clipping

![Soft Clipping](https://i.imgur.com/Bg2ygds.png "Soft Clipping")

Soft clipping, on the other hand, means that the peaks of a waveform experience a reduced increase, resulting in a more "rounded" appearance. Soft clipping in an op-amp occurs when the output signal approaches the supply voltage limits but does not abruptly cut off; instead, it gradually compresses the signal as it nears those limits. This typically happens when the input signal exceeds the linear operating range of the op-amp, causing the output to start saturating.

Harmonic Distortion in Operational Amplifiers

The Angstrom Method determines thermal conductivity and diffusivity by sending temperature pulses through a material and measuring at two points.

**The Angstrom Method is a method to determine the thermal conductivity or the thermal diffusivity of different materials. The method builds upon sending temperature pulses from a reservoir and measuring the temperature at different points in the material.**

## The Heat Equation

The situation is determined by the heat equation,

$$
    \frac{\partial T}{\partial t}(\mathbf{r},t) = \alpha \Delta T(\mathbf{r},t) - \epsilon T(\mathbf{r},t),
$$

where $T(\mathbf{r},t)$ is the temperature as a function of position $\mathbf{r}$ and time $t$, $\alpha$ is the thermal diffusivity of the material (<a href="https://www.amazon.com/Concepts-Thermal-Physics-Stephen-Blundell/dp/0199562105?crid=3M2TACZV99FLQ&dib=eyJ2IjoiMSJ9.1Upe8Fxd9wr9d7I6Ti7RFgLIEhwPL8AaTT0sSahT5HAD8Riv3p-qDFkupz2j8Mg6KQZIZ9HkaBy2zmTvbPcqpythdKGXNEUMUVUPwfybEL6mPFnbyquUB-ME3Zc-1ZSrms_TfZ2yLUheByWpQwt9YJoKvTVthEmYGwvwoGc9zNcyBKMojMJW_uSniMs5jHudcd6okW5zYHSrZ08WnqNdNzySm8OZ5oKTOLRSthLTQ6I.RCCeTr6wZbRMNcUVN9-sU0Alnot7tdHPY3dSjh78PYc&dib_tag=se&keywords=concepts+in+thermal+physics&qid=1733580143&sprefix=angstrom+methodconcepts+in+thermal+physics%2Caps%2C389&sr=8-1&linkCode=ll1&tag=engui-20&linkId=24019e7d94d809c8f3a0e8ed9e59d5f8&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 89</a>) and $-\epsilon T$ represents the heat losses to the ambient air. In the case of the material being an oblong rod, the situation can be approximated as one-dimensional, which gives the following formulation of the heat equation

$$
    \frac{\partial T}{\partial t}(x,t) = \alpha \frac{\partial^2 T}{\partial x^2}(x,t) - \epsilon T(x,t).
$$

Furthermore, $\alpha = \frac{\kappa}{c \rho}$, where $\kappa$ is the thermal conductivity, $c$ is the speific heat capacity and $\rho$ is the density, as well as $\epsilon = \frac{R O}{S c \rho}$ where $R$ is an emission constant, $O$ is the circumference of the rod and $S$ the cross-sectional area of the rod, meaning

$$
    \frac{\partial T}{\partial t}(x,t) = \frac{\kappa}{c \rho} \frac{\partial^2 T}{\partial x^2}(x,t) - \frac{R O}{S c \rho} T(x,t)
$$

According to the Angstrom method, a heater shall now be connected to one of the rod's ends, say at $x=0$, at which it works as a reservoir securing the boundary condition $T(0,t) = T_0 \cos \omega t$ for any chosen temperature $T_0$ och ''angular velocity'' $\omega$. Due to diffusion and emission, the temperature $T$ will not only vary sinusoidal with $t$ and $x$ but also exponentally decrease with $x$. Therefore, a solution with the form

$$
    T(x,t) = A e^{-ax} \cos(\omega t - b x),
$$

where $A$, $a$ (probably $>0$) and $b$ are constants, will be assumed. If this is put into the heat equation for the oblong rod then

$$
    -\omega A e^{-a x} \sin(\omega t - b x) =
$$
$$
 \frac{\kappa}{c \rho} A e^{-a x} \big((a^2-b^2) \cos(\omega t - b x) - 2 a b \sin(\omega t - b x)\big) 
$$
$$
- \frac{R O}{S c \rho} A e^{-a x} \cos(\omega t - b x),
$$

is required, which if the exponential factor $A e^{-a x}$ is canceled, gives two equations in $\sin$ and in $\cos$

$$
    \sin \implies -\omega = 2 a b \frac{\kappa}{c \rho},
$$

$$
    \cos \implies \frac{\kappa}{c \rho} (a^2-b^2) - \frac{R O}{S c \rho} = 0.
$$

## Using Two Measuring Points to Determine Thermal Conductivity

Furthermore, the Angstrom method is based on measuring the temperature, $T_1$ and $T_2$, at two points, $x_1$ and $x_2$ (with $x_2>x_1$ for simplicity), and by comparing the amplitudes, say $A_1$ and $A_2$ respectively, at those points. With the nature of the temperature function, it follows that

$$
    \frac{T_2}{T_1} = \frac{A_2}{A_1} = e^{-a(x_2-x_1)},
$$

or 

$$
    a = \frac{\log(\frac{A_1}{A_2})}{x_2-x_1},
$$

and by comparing the phases, say $\varphi_1$ and $\varphi_2$ respectivaly, it follows that

$$
    \varphi_2-\varphi_1 = b(x_2-x_1),
$$

or

$$
    b = \frac{\varphi_2-\varphi_1}{x_2-x_1}.
$$

and with these equations, $a$ and $b$ can be eliminated, which gives the following expression for the thermal conductivity $\kappa$

$$
    \kappa = \frac{c \rho \omega (x_2-x_1)^2}{2 (\varphi_1-\varphi_2) \log(\frac{A_1}{A_2})}.
$$

## Determining Thermal Diffusivity

Using the above equation for $\kappa$ we can insert that into the equation for the thermal diffusivity

$$
    \alpha = \frac{\kappa}{c \rho} = \frac{\omega (x_2-x_1)^2}{2 (\varphi_1-\varphi_2) \log(\frac{A_1}{A_2})}.
$$

So that is how you determine the thermal diffusivity and thermal conductivity of a material using the Angstrom Method.

Determining Thermal Conductivity Using The Angstrom Method

To apply models like the ARMA model to time series sample data, the seasonality and trend have to be removed. Let's look at how to do that.

**To apply models like the ARMA model to time series sample data, the seasonality and trend have to be removed. Let's look at how to do that.**

Let us say we have time series $X$ which has both a trend and a seasonality then we can do a decomposition of the time series and write it as:

$$
    X_t = m_t + s_t + Y_t,
$$

where $m_t$ is the trend component, $s_t$ is the seasonal component with period $d$ and $Y_t$ is the rest which is [weakly stationary](https://www.engineeringguidance.com/time-series-white-noise) and has mean of 0. Since we have singled out a seasonal component we know that $s_t = s_{t+d}$ because $d$ is the period of the seasonality.

Before we discuss how to eliminate the trend and seasonal component, we have to examine a couple of time series operators. 

## The Lag Operator

The lag operator $B$, also called the backshift operator, is defined as:

$$
    B X_t = X_{t-1}.
$$

So when taking the lag operator to the power of $k$ we get:

$$
    B^k X_t = X_{t-k}.
$$

## The Difference Operator

The difference operator $\nabla$ is defined as:

$$
    \nabla X_t = (1-B) X_t =  X_t - X_{t-1}.
$$

So when taking the difference operator to the power of $k$ and using the binomial theorem we get:

$$
    \nabla^k X_t = (1-B)^k X_t = \sum_{i=0}^k {k \choose i} (-1)^i B^i X_t 
$$
$$
= \sum_{i=0}^k {k \choose i} (-1)^i X_{t-i}.
$$

## The Lag-d Differencing Operator

The lag-d differencing operator $\nabla_d$ is defined as:

$$
    \nabla_d X_t = (1-B^d) X_t =  X_t - X_{t-d}.
$$

## Removing Seasonality in a Time Series

To eliminate seasonality in a time series we will have to use the lag-d differencing operator. Since the seasonality has a period of $d$ then:

$$
    \nabla_d X = m_t - m_{t-d} + s_t - s_{t-d} + Y_t - Y_{t-d} 
$$
$$
= \nabla_d m_t + \nabla_d Y_t,
$$

because $s_t = s_{t-d}$.

## Removing Trend in a Time Series

Eliminating trend in a time series is a bit more tricky than removing seasonality. We have to assume that

$$
    m_t = \sum_{i=0}^q a_i t^i,
$$

where $a_i$ is some real constant and $q < n$ where $n$ is the sample size of the data. Lets apply the difference operator to the power of $q$ on $m_t$:

$$
    \nabla^q m_t = \sum_{i=0}^q {q \choose i} (-1)^i \sum_{j=0}^q a_j (t-i)^j.
$$

Since the sums are finite we can change the order of them, so we get the following

$$
    \sum_{j=0}^q a_j \sum_{i=0}^q {q \choose i} (-1)^i (t-i)^j.
$$

Now, we focus on the inner sum:

$$
    \sum_{i=0}^q {q \choose i} (-1)^i (t-i)^j.
$$

By applying the binomial theorem and properties of polynomial expansions on $(t-i)^j$ we get:

$$
    (t-i)^j = \sum_{k=0}^j {j \choose k} t^{j-k} (-i)^k.
$$

By substituting this into our sum and changing order of the sums we get:

$$
    \sum_{i=0}^q {q \choose i} (-1)^i (t-i)^j 
$$
$$
= \sum_{k=0}^j {j \choose k} t^{j-k} (-1)^k \sum_{i=0}^q {q \choose i} (-1)^i i^k.
$$

Using the fact that

$$
    \sum_{i=0}^q {q \choose i} (-1)^i i^k = 0 \text{ for } k < q
$$

we find that the entire expression evaluates to 0 when $j < q$. When $j = q$ this inner sum evaluates to $q!$. Therefore, the final expression simplifies to:

$$
    a_q q!
$$

and this means that

$$
    \nabla^q m_t = a_q q!.
$$

So when we apply $\nabla^q$ to $X_t$ we get

$$
    \nabla^q X_t = a_q q! + \nabla^q Y_t.
$$

Given that $Y$ is regarded as a stationary process with a mean of zero, it can be demonstrated that this property extends to $\nabla^q Y$. Consequently, if $s = 0$, $\nabla^q X$ behaves as a stationary process with a mean of $a_q q!$.

Time Series Seasonality and Trend

The patch antenna connects the TRM's amplifiers to free space, transmitting electromagnetic waves with minimal loss and directionality for longer range.

**The patch antenna is the component used to connect the [Transmit Receive Module's](https://www.engineeringguidance.com//what-is-a-transmit-receive-module) amplifiers to free space, allowing for electromagnetic waves to be transferred with minimal losses. For most antennas to work efficiently, however, the radiated electromagnetic energy not only needs to be transmitted but also directed in order to be able to transmit efficiently and over longer distances.**

## What is an Antenna Array?

A single antenna element usually has a wide radiation pattern but for many applications such as for radars, it is necessary to have antennas with a high directivity that can create beams that focus the radiated power. To do this the most common way is to form an array of antenna elements, known as an antenna array. The total radiated field of such an array is determined by how the individual elements' radiation interfere in different directions which results in constructive interference in some directions while others experience destructive interference. 

## Two Element Antenna Array

![Two element linear array](https://i.imgur.com/CruAuUo.png "Two element linear array")

The smallest antenna array possible is a simple two element linear array as the one shown in the figure above consisting of two infinitesimal horizontal dipoles positioned at $\pm d/2$ along the z-axis. Assuming no coupling between the elements the radiated field is given by adding the field of the two individual antenna elements giving

$$
    \mathbf{E}_t = \mathbf{E_1} + \mathbf{E_2} 
$$
$$
= \hat{\theta}j\eta\frac{kI_0l}{4\pi}\left( \frac{e^{-j[kr_1-(\beta/2)]}}{r_1}\cos{\theta_1} + \frac{e^{-j[kr_2+(\beta/2)]}}{r_2}\cos{\theta_2}\right).
$$

where $\eta$ is the free space impedance, $k$ is the wavenumber, $I_0$ is the current amplitude, $l$ is the length of the dipoles, and $\beta$ is the phase of the signal. With the assumptions that $\theta\simeq\theta_1\simeq\theta_2$, the phase variations $r_1\simeq r-\frac{d}{2}\cos{\theta}$, $r_2\simeq r+\frac{d}{2}\cos{\theta}$ and the amplitude variation  $r_1\simeq r_2 \simeq r$ the equation reduces to 

$$
    \mathbf{E_t} = \hat{\theta}j\eta\frac{kI_0le^{-jkr}}{4\pi r}\cos{\theta}\cdot\left\{ 2\cos\left[\frac{1}{2}(kd\cos{\theta}+\beta)\right] \right\}.
$$

The last factor within curly brackets is called the array factor (AF). It can be seen by comparing the first and second equations that the radiated field from the array is the product of the radiated field from a single element multiplied by the AF. By varying the separation $d$ and/or the phase $\beta$ the total field of the array can be controlled. Most antennas however have a fixed separation by design and what is left to control is the phase. When controlling the phase for the antenna array it is then called a phased antenna array. 


![Beam angle antenna](https://i.imgur.com/zw4w0mC.png "Beam angle antenna")

The same principle as for the two element array still applies for arrays when extending it to higher dimensions and adding more elements. Using more elements can get even better directivity allowing for beams that can reach longer. An illustration of the beam angles $(\theta,\varphi)$ in three dimensions can be seen in the figure above where it can be seen that at $\theta = 0$ and $\varphi = 0$ the antenna array is sending radiation straight ahead. To get more accurate results, however, the interaction between the elements, known as mutual coupling, has to be accounted for.

## Mutual Coupling in Antenna Array

Mutual coupling in an antenna array is a term that is used to describe the interaction between antenna elements. It occurs mainly as a consequence of nearby elements receiving some of the radiated energy but can also occur because of reflections, reradiation, and scattering of electromagnetic radiation. This phenomenon exists for all types of arrays but becomes even more important to account for when working with large arrays. What happens is that when an element is excited by a current, the surrounding antenna elements will receive some of the energy leading to changes in their own respective excitation. This sums up to a total change in the excitation for all elements which will change the characteristics of the antenna seen from the driving source. One key factor is that the antenna's terminal impedance, as seen by the source, changes which might cause an impedance mismatch resulting in a decrease in radiated energy. While designing an antenna array the coupling is accounted for and the antenna array is typically impedance matched for a beam angle straight ahead. To assess the beamforming capability, however, all angles of interest have to be evaluated since the coupling tends to be worse for large beam angles.

## The Patch Antenna

A patch antenna is a type of microstrip antenna that is lightweight, easy to manufacture, and commonly used for a variety of applications. It typically consists of three layers; a ground plane, a substrate, and the patch elements. A patch antenna designed on a printed circuit board (PCB) requires a large grounded plane below the antenna elements with the substrate in between. The signal can be fed to the elements in different ways such as through a microstrip feed line or with a coaxial cable with its feed pin going through the substrate and soldered to the patch, the latter being shown in the figure above.

When designing a patch antenna array, lots of different aspects have to be taken into account. For example, the shape of the patch, the type of material, and the dimensions all have significant effects on the antenna array performance. 

## Design of Patch Antenna Array

While the layout of the patch can take different shapes to better suit a specific application, a rectangular patch antenna with probe feeding is shown in the figure above because of its simplicity. This shape has also proven to have better overall characteristics compared to other simple shapes such as triangular and circular patches.

For a large antenna array, most of the antenna elements will behave similarly with the exception of the elements close to the array edge. Therefore a practical way to design an antenna array is to use a unit cell representing one element within an infinite array. Using a unit cell reduces the simulation time and computational complexity significantly. It is a decomposition-based finite element technique for modeling repetitive antenna arrays and can be used because of the similarities of the electromagnetic fields for elements inside a large array. 

![Hex cell patch antenna](https://i.imgur.com/hgFSM7X.png "Hex cell patch antenna")

The figure above shows a hexagonal unit cell of one of the patch antennas in an infinite array. A hexagonal unit cell can fit more elements within the same amount of area compared to a square cell. The element spacing in all three directions is chosen to be half a wavelength ($\lambda/2$). Using $\lambda/2$ spacing of the elements is chosen since it gives a good trade off between directivity and stearability with low risks of grating lobes. Grating lobes are beams that can be created for larger scan angles if the spacing is higher which would cause a decrease in transmitted power in the desired direction.

Depending on the frequency at which the antenna is designed the remaining parameters such as cell width, patch width, feed pin offset, feed pin radius, and the feed pin's teflon thickness for the antenna element can optimized for this single frequency using a finite element analysis software and common measurements from manufacturers.

Patch Antenna Arrays

In fluid dynamics, there are many concepts, formulas, and coefficients. In this post the most fundamental ones are explained

**In fluid dynamics, there are many concepts, formulas, and coefficients. In this post the most fundamental ones are explained**

## Bernoulli's Principle

In incompressible isentropic flow, i.e. a flow where a negligible amount of energy is lost to non-reversible processes such as heat in turbulent flows, one can derive Bernoulli's principle from the conservation of energy

$$
    \frac{u^2}{2}+\frac{p}{\rho}+z\cdot g= \text{constant},
$$

where $z$ is the vertical position of a fluid element, this law states that an increase in velocity implies a decrease in pressure and vice versa. 

## Dimensionless parameters
Dimensionless parameters are used to make a parameter comparable with parameters of other types of flow. Below the relevant parameters are presented.

### Reynolds Number Formula
The Reynolds number, Re, is a dimensionless number that tells the relation between the inertia and viscosity in a newtonian fluid and is defined by the following equation

$$
    \text{Re} = \frac{\rho V L}{\mu} = \left\{ \frac{\text{Inertia}}{\text{Viscosity}} \right\}.
$$

$V$ and $L$ are the characteristic velocity and length of the flow. A high Re is associated with a fast, large-scale turbulent flow, whereas a low Re corresponds to a slow, viscous flow. Fast flows of gas imply a relatively high Re.

### Drag Coefficient Formula

The drag coefficient of a body, $C_D$, is a dimensionless number indicative of drag force. $C_D$ is defined by the following equation

$$
    C_D = \frac{F_D}{\frac{1}{2}\rho u_{\infty}^2 A},
$$

where $A$ is the projected area from the front of the body and the denominator $(\frac{1}{2} \rho  u_{\infty}^{2})$ is the dynamic pressure of the free stream.

### Lift Coefficient Formula

The lift coefficient of a body, $C_L$ is defined by the following equation

$$
    C_L = \frac{F_L}{\frac{1}{2}\rho u_{\infty}^2 A},
$$

where $A$ is the projected area from the front of the body. $C_L$ is a dimensionless coefficient that increases with lift force. Depending on the area of study, $C_L$ will have different directions. For example, when designing race cars, $C_L$ will indicate the lift in the negative $ z$ direction, i.e., a high downforce will correspond to a large $C_L$. When designing airplanes it will be in the other direction to indicate lift.

### Pressure Coefficient Formula

The pressure coefficient, $C_p$, is a dimensionless parameter that describes the relative pressure. $C_p$ is defined as 

$$
    C_{p}=\frac{p-p_{\infty}}{\frac{1}{2} \rho  u_{\infty}^{2}},
$$

where $p$ is the static pressure at the point where $C_p$ is calculated, while $p_{\infty}$ is the static pressure in the free stream. 

### Skin Friction Coefficient Formula

The skin friction coefficient, $C_f$, is a dimensionless parameter defined as

$$
    C_{f}=\frac{\tau _{w}}{\frac{1}{2} \rho  u_{\infty}^{2}}.
$$

### Vorticity Formula

The vorticity, $\vec{\zeta}$, is defined as the curl of the velocity field,  $\nabla\,\times\,\vec{u}$, which relates to the rotation of a velocity field. If the vorticity is zero the flow has no rotation and is called irrotational. The vorticity in the $x$-direction is given by 

$$
    \zeta_x= \frac{\partial w}{\partial y} - \frac{\partial v}{\partial z}.
$$

## Navier-Stokes Equations

A fluid can be described with the Navier-Stokes equations if it is incompressible and newtonian as 

$$
    \rho \left( \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla) \vec{u} \right) = \rho \vec{g} - \nabla p + \mu \Delta \vec{u}. 
$$

### RANS Formula

In turbulent flow, some fluctuations cause rapid changes in the velocity and pressure in the Navier-stokes equations. These rapid changes can be easier to account for if the velocity and pressure are split into a mean variable $\bar{u}$ and a fluctuation variable $u'$. The velocity in the $x$-direction is therefore $u = \bar{u} + u'$. If those splits are put into the Navier-stokes equations we get the Reynold's Averaged Navier Stokes (RANS) equations after some derivations. The RANS in $x$-direction is 

$$
     \rho\frac{d \overline{u}}{dt}=-\frac{\partial \overline{p}}{\partial x}+\rho g_x+\frac{\partial }{\partial x}\left(\mu\frac{\partial \overline{u}}{\partial x}-{\rho \overline{u^{\prime 2}}}\right)
$$
$$
+\frac{\partial }{\partial y}\left(\mu\frac{\partial \overline{u}}{\partial y}-{\rho \overline{u^{\prime}v^{\prime}}}\right)+\frac{\partial }{\partial z}\left(\mu\frac{\partial \overline{u}}{\partial z}-{\rho \overline{u^{\prime}w^{\prime}}}\right).
$$
 
## Flow Past Boundary

Bodies, such as plates, immersed in a fluid stream create a boundary layer flow. The boundary layer is defined as the region where the flow velocity is less than 99 percent of the velocity of the external flow ($u_\infty$). The boundary layer flow is made up of two parts, the laminar part and the turbulent part:

$$
    \frac{\delta}{x} \approx
    \left\{
    \begin{aligned}
        & \frac{5}{\text{Re}_x^{1/2}} \qquad 10^3 < \text{Re}_x < 10^6 \text{ laminar flow} \\
        & \frac{0.16}{\text{Re}_x^{1/7}} \qquad 10^6 < \text{Re}_x \text{ turbulent flow}
    \end{aligned}
    \right. .
$$

$\delta$ is the thickness of the boundary layer and $\text{Re}_x$ the Reynolds number at $x$. An illustration of he boundary level flow can be seen in figure below.


![Boundary layer flow](https://i.imgur.com/HpCZOT0.png "Boundary layer flow")

### Flow Separation

Backflow, i.e. flow going in the opposite direction of the freestream, at the wall indicates that the flow is separated. The point of separation is when the wall shear is 0 which means the velocity gradient is 0.

### The Logarithmic Overlap Law

The turbulent flow near a wall can be broken up into three regions based on which shear stress dominates:

* Wall layer: Viscous shear dominates
* Outer layer: Turbulent shear dominates
* Overlap layer: A mix of the shear types

The inner layers velocity can be modeled as $u = f(\mu, \tau_w, \rho, y)$, $\tau_w$ is the wall shear stress. By dimensional analysis the following relation can be deducted:

$$
    u^+ = \frac{u}{u^*} = F \left(\frac{y u^*}{\nu} \right),
$$

where $u^* = \sqrt{\tau_w/\rho}$. An akin method can be done with the outer layer, $(u_\infty - u)_{\text{outer}} = g(\delta, \tau_w, \rho, y)$ and the dimensional analysis leads to:

$$
    \frac{u_\infty - u}{u^*} = G\left( \frac{y}{\delta} \right).
$$

The overlap layer makes it so these two layers overlap smoothly and for this to be true the overlap velocity must vary logarithmically with y:

$$
    \frac{u}{u^*} = \frac{1}{\kappa} \ln{\frac{y u^*}{\nu}} + B,
$$

where $\kappa \approx 0.41$ and $B \approx 5.0$. 
We define the dimensionless number $y^+$ as

$$
    y^+ =\frac{y u^*}{\nu},
$$

where $u^* = \sqrt{\tau_w/\rho}$. For $y^+ \in [0,5]$ the inner layer dominates and is proportional to $y^+$. It can be described with the dimensionless variable

$$
    u^+ = \frac{u}{u^*} = y^+.
$$

For $y^+ \in [30, 10^3]$ the overlap layers dominates and can be described as

$$
    u^+ = \frac{1}{\kappa} \ln{y^+} + B,
$$

where $\kappa \approx 0.41$ and $B \approx 5.0$. For $y^+ \in [5,30]$ the inner layer has to curve so it merges with the overlap layer. Only experimental data exists for this interval making it difficult to model.

Fundamental Fluid Dynamics Theory

A varactor diode is a diode with adjustable capacitance controlled by the reverse voltage, functioning as a variable capacitor.

**What is a varactor diode? The varactor diode, which stands for variable capacitor, is a diode with a well defined capacitance that can be controlled with the reverse voltage applied to it.**

## How Does a Varactor Diode Work?

The varactor's depletion region width $d$, seen in the illustration of the varactor diode in the figure above, increases when the reverse voltage over the diode increases and can be likened to the distance between the plates of a plate capacitor being increased as seen in 

$$
    C = \frac{A \varepsilon}{d},
$$

where $C$ is the capacitance, $A$ is the cross section area of the plate capacitor, $\varepsilon$ the dielectric constant and $d$ is the distance between the two plates. The increase of $d$ due to the increase of reverse voltage causes the capacitance of the varactor to decrease. These characteristics enable the varactor's capacitance to be controlled and set to a specific value. All varactor diodes do however have a maximum allowed reverse voltage known as the breakdown voltage. Exceeding the breakdown voltage could break the varactor which is why the reverse voltage must be kept below this limit.

## Varactor Diode Symbol

![Varactor Diode Symbol](https://i.imgur.com/sTzVizV.png "Varactor Diode Symbol")

The varactor's circuit can be seen in the figure above. As can be seen, it is a combination of the symbol of a diode and a capacitor.

## Varactor Diode Equivalent Circuit

The figure below shows the SPICE (Simulation Program with Integrated Circuit Emphasis) model used for modeling the varactor as an equivalent circuit. It is derived from *Skyworks Solutions Inc*'s datasheets. The SPICE model contains parasitic components $C_p$ and $R_s$ and the package inductance $L_s$. The difference compared to *Skyworks Solutions Inc*'s datasheet is the change from the diode model to a simple variable capacitance $C_\text{var}$.

![Linear SPICE model of varactor diode](https://i.imgur.com/Yj6lro4.png "Linear SPICE model of varactor diode")

## Nonlinear Varactor Diode Modeling

While the SPICE model shown in the figure above works for simulating S-parameters, a nonlinear model that accounts for the voltage swing across the varactor is required if more accurate results are wished for. The circuit simulation software *Advanced Design System* provides a solution for this using a built-in component called Symbolically Defined Devices (SDD). The SDD enables the creation of equation-based nonlinear components that are a function of the port voltages, currents, and their derivates. In the figure below an alternative schematic using the SDD can be seen.

![Nonlinear SPICE model of varactor diode](https://i.imgur.com/549EoOt.png "Nonlinear SPICE model of varactor diode")

To model the nonlinear varactor using the SDD, an equation describing the capacitance as a function of reverse voltage $C(v)$ is required. This can either be obtained by fitting a curve to measured data or, if available, be found in the datasheet of the varactor. Integrating $C(v)$ with respect to $v$ gives the charge $Q(v)$ as a function of reverse voltage which can be used as input for the SDD.

The Varactor Diode - Symbol and Equivalent Circuit

A Transmit Receive Module handles transmission and reception functions in a phased array before the signal is emitted through the antenna element.

**A Transmit Receive Module carries out the necessary functions in a phased array just prior to the signal being emitted through the antenna element. It encompasses both transmission and reception capabilities, which is why it’s called a Transmit Receive Module.**

A Transmit Receive Module's (TRM) main components are amplifiers, filters, phase shifters, switches, and a circulator. A schematic for a typical TRM is shown in the figure. The TRM can be switched into either be in transmit or receive mode which alters the path of the signal. Some of the functionality of the TRM is used both while transmitting and receiving, namely the phase shifter and attenuator. The phase shifter is what enables the beamforming to scan the beam in different directions while the attenuator is used to control the signal amplitude. Sharing these components between both modes has several benefits. Not only does it reduce the number of active components that would increase the complexity but it also saves weight and volume. While the savings for a single TRM might be low, the total decrease in weight and volume for an array of thousands of TRMs is significant.

## Power Amplifiers in TRMs
Inside the TRM there are two main amplifiers which are shown in the figure above. In the receive path sits the low noise amplifier (LNA) that is in charge of amplifying the received signal from the antenna. This very weak signal is greatly affected by noise which is why minimizing the noise that the amplifier adds is of high importance. In addition to the LNA, a receive protection circuit is used to make sure that signals with too high power do not get into the LNA. 

In the transmit path sits the HPA. It is responsible for amplifying the signal to the desired power level and ideally provides perfectly linear amplification. However, to maximize the transmitted power a trade off between linearity and output power usually has to be made. For example, having too much output power will start to compress the output signal causing distortions and these are things an HPA designer has to take into account. Furthermore, an HPA is designed to be presented to a specific impedance, both on the input and output ports. Typically both these are set to 50 $\Omega$ by design, requiring the rest of the system to be designed around the same impedance but this varies for different applications.

In the case of an HPA with a 50 $\Omega$ output impedance that is presented to a load differing from the designed one, reflections will occur. These reflections can greatly affect the HPA resulting in reduced power efficiency and gain which can have a severe impact on its performance or potentially even break it. Making sure that the HPA is protected from reflections is therefore of high importance to achieve greater and more stable performance.

## The TRM's RF Circulator

Between the TRM and antenna, a circulator is most commonly used. A circulator is a passive component and a type of duplexer, allowing for bidirectional communication of a single path. In practice this means that signals coming from the HPA in the transmit path will be sent to the antenna and signals coming from the antenna will be sent to the receive path, thereby separating the signals without using a switch to select one path. Additionally, it also makes sure that reflections occurring at the antenna do not travel back to the HPA, degrading its performance. While the use of a circulator handles the main problems of connecting the TRM and antenna it also comes with the drawbacks of being expensive and requiring a considerable amount of volume. Another drawback is that while it typically adds less than 1 dB loss of signal power, the total loss of power can be larger when accounting for reflections between the circulator and antenna. 

## Transmit Receive Modules in Phased Antenna Arrays

The enabling technology behind Active Electronically Scanned Arrays (AESAs) comes from the ability to control hundreds or thousands of TRMs individually. Each TRM connects to an individual antenna element and is collectively controlled by a beamformer that decides how every TRM should behave at all times. This allows an AESA to create multiple directed beams which can be used for multiple tasks at once or to improve one further. 

The most common implementation of AESAs is for radar systems. Radar systems based on AESAs have several advantages compared to other technologies such as mechanically scanned arrays (MSAs) and passive electronically scanned arrays (PESAs). Both AESAs and PESAs have quicker scanning speeds than MSAs since they do not require any mechanically moving parts. 

Furthermore, the main difference between PESAs and AESAs is the number of TRMs. While a PESA has individual phase shifters, it only has a single TRM compared to AESAs which have one for each element. This gives many advantages such as allowing for the multibeam functionality of AESAs that allow radar systems to track a larger amount of targets. The individual beams can also be used to transmit multiple frequencies at once to disguise their signals which is a common task within electronic warfare and is not possible with PESAs that only can transmit a single beam. Another key feature of having individual TRMs is that even if a few TRMs fail, the functionality of the array remains which adds redundancy to the system. 

## Beamformer
Connecting to all the TRMs is the beamformer. The beamformer is what distributes and combines the signal between the TRMs in addition to steering the TRM's respective phase shift and amplitude at all times. When in transmission mode, the beamformer takes a modulated waveform and distributes it between the TRMs, and decides which phase shift should be added to each of them along with the amplitude to get the desired beam. In receive mode, the beamformer still controls the TRMs but instead of distributing a waveform, it is responsible for combining the received TRM signals into a single signal. Both in transmit and receive mode, the beamformer's phase control of the TRMs is crucial to get a beam that transmits or receives in the desired direction.

What is a Transmit Receive Module?

This guide covers identifying white noise and iid noise in residuals with ACF/PACF plots and the Ljung-Box test.

**In time series sample data there is always noise. But there exist different types of noise. In this post we will run through regular white noise and iid noise and how to identify them in your residuals with ACF/PACF plots and the Ljung-Box test after for example applying an [ARMA](https://www.engineeringguidance.com/the-arma-process) or GARCH model. First we need to define weak stationarity before we define white noise.**

## What is Weak Stationarity?

A time series $X$ is weakly stationary if

- $\mathbb{E}[X_t] = \mu = \text{const}$ for all $t \in \mathbb{Z}$
- $\text{Var}[X_t] = \sigma^2 = \text{const} < \infty$ for all $t \in \mathbb{Z}$
- $\text{Cov}[X_r, X_s] = \text{Cov}[X_{r+h}, X_{s+h}]$ for all $s,r,h \in \mathbb{Z}$

 (<a href="https://www.amazon.com/Introduction-Forecasting-Springer-Texts-Statistics/dp/3319298526?crid=3J0ML8ICAR1T&dib=eyJ2IjoiMSJ9.YtB4SEwuz6VeggvSidRxhqyID70S1JcDM9g5RnheKrzg-X_W3cD9LPYBYEKZ5XXE2ELLEd-goyZauiisohEzWeEqinBEFIqg7CcSVYAAKcHAXASTIyzmEFRkPMkQtpO4Ddfw0vlJHqvxewxJStEyG34itUpZXzgZ0708ELIbYXREcl1PSyHuzthWI49xI-8XIAkdPfSm-Xo2CAiVb2GjM2er6B_qvIUPWM-uGYCnCDA.VXMBJweWy62TXA3CkV-IZfnFaoHog4kziLyP-rHTnW4&dib_tag=se&keywords=introduction+to+time+series+analysis+and+forecasting&qid=1734113696&sprefix=introduction+to+time+series%2Caps%2C241&sr=8-5&linkCode=ll1&tag=engui-20&linkId=73e37b590b35d5100ac24e7e309c4a83&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 13</a>).

## Defintion of White Noise

A time series $X$ is said to be white noise if

- $X$ is weakly stationary
- $\gamma(h) = \text{Cov}[X_t, X_{t+h}] = 0$ where $h \neq 0$

(<a href="https://www.amazon.com/Introduction-Forecasting-Springer-Texts-Statistics/dp/3319298526?crid=3J0ML8ICAR1T&dib=eyJ2IjoiMSJ9.YtB4SEwuz6VeggvSidRxhqyID70S1JcDM9g5RnheKrzg-X_W3cD9LPYBYEKZ5XXE2ELLEd-goyZauiisohEzWeEqinBEFIqg7CcSVYAAKcHAXASTIyzmEFRkPMkQtpO4Ddfw0vlJHqvxewxJStEyG34itUpZXzgZ0708ELIbYXREcl1PSyHuzthWI49xI-8XIAkdPfSm-Xo2CAiVb2GjM2er6B_qvIUPWM-uGYCnCDA.VXMBJweWy62TXA3CkV-IZfnFaoHog4kziLyP-rHTnW4&dib_tag=se&keywords=introduction+to+time+series+analysis+and+forecasting&qid=1734113696&sprefix=introduction+to+time+series%2Caps%2C241&sr=8-5&linkCode=ll1&tag=engui-20&linkId=73e37b590b35d5100ac24e7e309c4a83&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 14</a>).

## Definition of IID Noise

What is iid noise? A time series $X$ is said to be iid noise if it meets the conditions for regular white noise above but also $X_r$ and $X_s$ ($r \neq s$) is independent of each other and identically distributed for all $r,s \in \mathbb{Z}$, that is why it is called iid noise where iid stands for independent and identically distributed. 

Before going through the tricks to identifying the two different kind of noises we need to go through the definitions of the autocorrelation function (ACF) and the partial autocorrelation function (PACF).

## The Autocorrelation Function (ACF)

The autocorrelation function for the time series $X$ is defined as:

$$
    \rho_X(h) = \frac{\gamma(h)}{\gamma(0)} = \frac{\text{Cov}[X_t, X_{t+h}]}{\text{Var}[X_t]}.
$$

## The Partial Autocorrelation Function (PACF)

The autocorrelation function $\alpha(h)$ for the time series $X$ is defined by $\alpha(0) = 1$ and $\alpha(h) = \phi_{hh}$ for $h \geq 1$ where $\phi_{hh}$ is the last component of

$$
    \phi_h = ([\gamma(i-j)]_{i,j=1}^{h})^{-1} [\gamma(1), ... , \gamma(h)]'.
$$

## ACF Plot Interpretation - How to Identify White Noise

The first simple thing you could do to see if your data is just white noise is if it looks like it has no structure. The second thing is looking at the sample ACF and PACF graphs for the data. Lets look at the definition of white noise and the ACF. We know that $\gamma(h) = 0$ for $h \geq 1$ for white noise which means $\rho_X(h) = 0$ for $h \geq 1$. So all lags greater than zero should be within the red striped confidence band ($\pm 1.96/\sqrt{n}$, where $n$ is the sample size) of the ACF (<a href="https://www.amazon.com/Introduction-Forecasting-Springer-Texts-Statistics/dp/3319298526?crid=3J0ML8ICAR1T&dib=eyJ2IjoiMSJ9.YtB4SEwuz6VeggvSidRxhqyID70S1JcDM9g5RnheKrzg-X_W3cD9LPYBYEKZ5XXE2ELLEd-goyZauiisohEzWeEqinBEFIqg7CcSVYAAKcHAXASTIyzmEFRkPMkQtpO4Ddfw0vlJHqvxewxJStEyG34itUpZXzgZ0708ELIbYXREcl1PSyHuzthWI49xI-8XIAkdPfSm-Xo2CAiVb2GjM2er6B_qvIUPWM-uGYCnCDA.VXMBJweWy62TXA3CkV-IZfnFaoHog4kziLyP-rHTnW4&dib_tag=se&keywords=introduction+to+time+series+analysis+and+forecasting&qid=1734113696&sprefix=introduction+to+time+series%2Caps%2C241&sr=8-5&linkCode=ll1&tag=engui-20&linkId=73e37b590b35d5100ac24e7e309c4a83&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 16</a>) in the plot below. So the following sample ACF plot indicates white noise:

![Sample ACF of the sample data](https://i.imgur.com/3exyzTL.png "Sample ACF of the sample data")

The same can be said for the PACF. If $X$ is white noise then the vector $[\gamma(1), ... , \gamma(h)]'$ will just be filled with zeros and thus $\alpha(h) = 0$ for $h \geq 1$.

## How to Identify IID Noise

To identify iid noise we will have to look at the sample ACF and PACF of the data squared, $X^2$. Lets start with the ACF. We have:

$$
    \gamma_{X^2}(h) = \text{Cov}[X^2_t, X^2_{t+h}]
$$
$$
 = \mathbb{E}[(X_t^2 - \mathbb{E}[X_t^2])(X_{t+h}^2 - \mathbb{E}[X_{t+h}^2])]
$$

and now we use the property of iid noise that $X_t$ and $X_{t+h}$ is independent of each other so we get the following:

$$
    (\mathbb{E}[X_t^2 - \mathbb{E}[X_t^2]]) (\mathbb{E}[X_{t+h}^2 - \mathbb{E}[X_{t+h}^2]])
$$
$$
 = (\mathbb{E}[X_t^2] - \mathbb{E}[\mathbb{E}[X_t^2]]) (\mathbb{E}[X_{t+h}^2] - \mathbb{E}[\mathbb{E}[X_{t+h}^2]]).
$$

Lastly we use the general property that $\mathbb{E}[\mathbb{E}[Y]] = \mathbb{E}[Y]$ so we get:

$$
    (\mathbb{E}[X_t^2] - \mathbb{E}[X_t^2]) (\mathbb{E}[X_{t+h}^2] - \mathbb{E}[X_{t+h}^2]) = 0.
$$

So the sample ACF for $X^2$ should be in the confidence band for all lags greater than zero. For the PACF we can see that the vector $[\gamma_{X^2}(1), ... , \gamma_{X^2}(h)]'$ would be zero for all $h$. Therefore the PACF should also be within the confidence band for lags greater than zero if $X$ were to be iid noise which we can see in the following plot:

![Sample PACF of the sample data squared](https://i.imgur.com/4ywfylN.png "Sample ACF of the sample data squared")

So to conclude, for iid noise $X$ and $X^2$'s sample ACF and PACF should all have to be within the confidence band for lags greater than zero.

## Ljung-Box Test for IID Noise

Lastly one can do a Ljung-Box test to determine whether the data is independently distributed. The test is defined by the following hypothesis:

$$
        \begin{aligned}
            & H_0 : X \sim \text{IID}(\mu, \sigma^2) \\
            & H_1 : X \not\sim \text{IID}(\mu, \sigma^2) 
        \end{aligned}
$$

The Ljung-Box test has a test statistic which is:

$$
    \lambda = n (n + 2) \sum_{i=1}^h \frac{\hat{\rho}(i)^2}{n-i},
$$

where $\hat{\rho}(i)$ is the sample ACF at lag $i$, $n$ is the sample size and $h$ is the number of lags being tested. The null-hypothesis, $H_0$, is rejected at the $\alpha$-level, $\alpha \in (0,1)$, if $\lambda > \chi_{1-\alpha, h}^2$, where $\chi_{1-\alpha, h}^2$ denotes the $1 - \alpha$-quantile of the $\chi^2$-distribution with $h$ degrees of freedom.

ACF and PACF of White Noise and IID Noise

The continuity equation in electromagnetism gives us a formula that relates the electrodynamics of the charge density to the current density in a volume.

**The continuity equation is an important tool that one can use whether it is in a [fluid dynamics](https://www.engineeringguidance.com/continuity-equation-for-fluids-in-different-coordinate-systems) problem or an electromagnetic problem. The continuity equation in electrodynamics relates the charge density to the current density in a volume.**

## Derivation of Continuity Equation For Charge and Current Density

Lets derive the continuity equation for electric charges. If you are just searching for it, this is what it looks like:

$$
    \frac{\partial \rho}{\partial t} = - \nabla \cdot \vec{j},
$$

where $\rho$ is the charge density and $\vec{j}$ is the current density (current divided by an area). Lets get a few things straight first. When we say the charge is conserved we are talking about the charge in an arbitrary volume. Additional charge can not be created from some sort of source inside the volume. The only way to change the charge is a flow of new charges from outside the volume, that would be our $\vec{j}$. So lets observe an infinitesimal surface element $d \vec{S}$ on the surface of our volume. The surface element is vector that is normal to the surface. So the infinitesimal charge $d q$ that passes through the surface element can be written as:

$$
    dq = \vec{j} \cdot d \vec{S} dt.
$$

We'll use this expression later. The total charge $q$ for the whole volume $V$ can be written as:

$$
    q = \int_V \rho dV,
$$

where earlier mentioned $\rho$ is the charge density for the volume. We can now use the conservation of charge and look at the change of the charge through time and flux of charge through the volume's surface $\partial V$. So we have two expressions, firstly:

$$
    \frac{dq}{dt} = \int_V \frac{\partial \rho}{\partial t} dV,
$$

which is just the derivative of the total charge for the volume. Secondly, we have:

$$
    \frac{dq}{dt} = - \int_{\partial V} \vec{j} \cdot d \vec{S},
$$

which uses the expression for the charge through a surface element. The reason for the minus sign is because we are calculating the inflow of charges but the surface element vector is pointing outwards from the volume and we have to counteract that. We are almost at our sought after expression but we'll have to use one more trick, the divergence theorem aka Guass's theorem which states:

$$
    \int_V \nabla \cdot \vec{F} dV = \int_{\partial V} \vec{F} \cdot d \vec{S},
$$

which fits perfectly with our charge flux expression:

$$
     - \int_{\partial V} \vec{j} \cdot d \vec{S} = - \int_V \nabla \cdot \vec{j} dV.
$$

So now we have:

$$
    \frac{dq}{dt} = \int_V \frac{\partial \rho}{\partial t} dV = - \int_V \nabla \cdot \vec{j} dV.
$$

As you can see we have two expression with a volume integral for the volume $V$. So we can remove the integral and we finally get our continuity equation for electric charges:

$$
    \frac{\partial \rho}{\partial t} = - \nabla \cdot \vec{j}
$$

(<a href="https://www.amazon.com/Introduction-Electrodynamics-David-J-Griffiths/dp/1009397753?crid=QV9D908FRQTE&dib=eyJ2IjoiMSJ9.Mb6GSnIftSBmmzB4_-XMmGFrU2c96t3YWJ8nI2Kjqn4Wvagg4awX0SQxHvVhGN7R9ncNrO9ZjVAzQK9VoJJxproGHnOTrQg0winQDPsmMJPc2mgqE71_M5Kh_dX4nctptwSmlc2w-NWOEdKOaItwdCOYgiTY_NCb5BAu0tlMIcxpehKISEgNMmrrrbknn3QqQ3UeP3Nq7irvPOZAPhIZ2UFewmWK0CXog7TOHkYHb8H-6QU7SrVoz1k7RcUC2oJT06XqnxZo8ZVfLFmss4AYlqBtNjefCr3JrtUUJ9vJsRI.aUhFL3L0ym9fRwv4CxGw7nOLwtwhWE8Igmx0DLoAVbA&dib_tag=se&keywords=electromagnetism&qid=1733583598&sprefix=electromagnetism%2Caps%2C208&sr=8-5&linkCode=ll1&tag=engui-20&linkId=a3d8de438cf500d32cb76df6633e921b&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 345</a>). Now lets see how the continuity equation looks in different coordinate systems.

## Continuity Equation in Cartesian Coordinates

The continuity equation for charges in Cartesian coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{\partial j_x}{\partial x} + \frac{\partial j_y}{\partial y} + \frac{\partial j_z}{\partial z} = 0
$$

(<a href="https://www.amazon.com/Handbook-Mathematics-I-N-Bronshtein/dp/3662462206?crid=2EY5M5H7FT2IH&dib=eyJ2IjoiMSJ9.DwfzsUXXkM7e94_ArqNkD7rI9LmiHe7Nl2zScHPEeCKOU1geuDYJusAIGYH94tYi0OQ9FNiQ4yi64-rEENwjhruC-ruHygf1MlPUoylnNrhbsTbUckL03fNH60wBfKaC-Uwew6hftZIJpPJZUSt1u1-MP39iufzwGsdI5WdbGMJm9l7J_Ea3MFirPBITTnUrW7Zx80gGJl0GXL1JdcZzsvBPWQ944WSWqiC9qxKd7fI.duuYoB-JYkro9TOd_-NwXkehbCumcwyC_VNTTpNfVbw&dib_tag=se&keywords=mathematics+handbook&qid=1733585380&sprefix=mathematics+hand%2Caps%2C193&sr=8-1&linkCode=ll1&tag=engui-20&linkId=4f7b8f54bfd4a3da9ba225ca7cf7a7f8&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 712</a>).

## Continuity Equation in Cylindrical Coordinates

The continuity equation for charges in cylindrical coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r}(r j_r) + \frac{1}{r} \frac{\partial}{\partial \theta}(j_\theta) + \frac{\partial}{\partial z}(j_z) = 0
$$

(<a href="https://www.amazon.com/Handbook-Mathematics-I-N-Bronshtein/dp/3662462206?crid=2EY5M5H7FT2IH&dib=eyJ2IjoiMSJ9.DwfzsUXXkM7e94_ArqNkD7rI9LmiHe7Nl2zScHPEeCKOU1geuDYJusAIGYH94tYi0OQ9FNiQ4yi64-rEENwjhruC-ruHygf1MlPUoylnNrhbsTbUckL03fNH60wBfKaC-Uwew6hftZIJpPJZUSt1u1-MP39iufzwGsdI5WdbGMJm9l7J_Ea3MFirPBITTnUrW7Zx80gGJl0GXL1JdcZzsvBPWQ944WSWqiC9qxKd7fI.duuYoB-JYkro9TOd_-NwXkehbCumcwyC_VNTTpNfVbw&dib_tag=se&keywords=mathematics+handbook&qid=1733585380&sprefix=mathematics+hand%2Caps%2C193&sr=8-1&linkCode=ll1&tag=engui-20&linkId=4f7b8f54bfd4a3da9ba225ca7cf7a7f8&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 712</a>).

## Continuity Equation in Spherical Coordinates

The continuity equation for charges in spherical coordinates is given by:

$$
    \frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 j_r) + \frac{1}{r \sin{\theta}} \frac{\partial}{\partial \theta}(j_\theta \sin{\theta}) + \frac{1}{r \sin{\theta}} \frac{\partial}{\partial \varphi}(j_\varphi) = 0
$$

(<a href="https://www.amazon.com/Handbook-Mathematics-I-N-Bronshtein/dp/3662462206?crid=2EY5M5H7FT2IH&dib=eyJ2IjoiMSJ9.DwfzsUXXkM7e94_ArqNkD7rI9LmiHe7Nl2zScHPEeCKOU1geuDYJusAIGYH94tYi0OQ9FNiQ4yi64-rEENwjhruC-ruHygf1MlPUoylnNrhbsTbUckL03fNH60wBfKaC-Uwew6hftZIJpPJZUSt1u1-MP39iufzwGsdI5WdbGMJm9l7J_Ea3MFirPBITTnUrW7Zx80gGJl0GXL1JdcZzsvBPWQ944WSWqiC9qxKd7fI.duuYoB-JYkro9TOd_-NwXkehbCumcwyC_VNTTpNfVbw&dib_tag=se&keywords=mathematics+handbook&qid=1733585380&sprefix=mathematics+hand%2Caps%2C193&sr=8-1&linkCode=ll1&tag=engui-20&linkId=4f7b8f54bfd4a3da9ba225ca7cf7a7f8&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 712</a>).

Continuity Equation in Electromagnetism For Current Density

The moment of inertia depends on the axis of rotation. Let's calculate the MOI formula for the middle and end of a rod and a homogenous and hollow disk.

**The moment of inertia depends on the axis of rotation. Let's calculate the MOI formula for the middle and end of a rod and a homogenous and hollow disk.**

We know that the moment of a point of inertia of a point mass is

$$
    M R^2,
$$

where $M$ is the mass and $R$ is the distance between the mass and the rotational axis. We can use this to calculate the moment of inertia for objects that is not point masses. We do this by integrating over all infinitesmal mass elements $dm$ of the object. So our moment of inertia is now

$$
    I = \int_0^M r^2 dm.
$$

## The Moment of Inertia of a Rod

The moment of inertia is dependent on the location of the rotational axis, therefore the moment of inertia will vary for a rod depending on the axis of rotation. We have our general formula for the moment of inertia

$$
    I = \int_0^M r^2 dm.
$$

We now have to identify what $dm$ is. Since the rod can be looked at in one dimension

$$
    dm = \frac{M}{L} dr
$$

where $L$ is the length of the rod. $\frac{M}{L}$ can be looked at as an one dimensional density for the rod.
### Moment of Inertia of a Rod About Its End
  Now lets calculate the moment of inertia if the rod is rotating at one of its ends. The moment of inertia then becomes

$$
    I = \int_0^L r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_0^L
$$
$$
 = \frac{M}{L} \frac{L^3}{3} = M \frac{L^2}{3}.
$$

### Moment of Inertia of a Rod About Its Center
 If the axis of rotation is in the middle of the rod the moment of inertia then becomes

$$
    I = \int_{-L/2}^{L/2} r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_{-L/2}^{L/2}
$$
$$
 = \frac{M}{3L} \left[ \frac{L^3}{8} - \left(\frac{-L^3}{8} \right) \right] = M \frac{L^2}{12}.
$$

### Moment of Inertia of a Rod General Formula
 Lets create a more general formula. Lets use the variable $x$ for the point of rotation. Our moment of inertia integral then becomes

$$
    I = \int_{-x}^{L-x} r^2 \frac{M}{L} dr = \frac{M}{L} \Bigg[ \frac{r^3}{3} \Bigg]_{-x}^{L-x}
$$
$$
 = \frac{M}{3L} \left[ (L-x)^3 + x^3 \right].
$$

## The Moment of Inertia of a Disk

![Image of the moment of inertia for disks](https://i.imgur.com/lyfEGNK.png "Moment of inertia for a homogeneous disk and a hollow disk")

Lets move on to calculate the moment of inertia for a disk (this also applies to a cylinder which is kind of an elongated disk).
### Moment of Inertia of a Disk Derivation
 The axis of rotation will be located in the center of the disk but we will get back to calculating the moment of inertia of a displaced rotational axis. Like usual we start of with the classic integral for the moment of inertia

$$
    I = \int_0^M r^2 dm.
$$

Defining $dm$ is a little bit more tricky than for the rod. In this problem we have to go two dimensional. Lets start with defining our two dimensional density which is

$$
    \frac{M}{\pi R^2}.
$$

The denominator is just the area of the disk. Our area element is 

$$
    r dr d\theta.
$$

I won't go into to much detail here, you'll have to do some repetition on polar coordinates if you don't understand how we get this area element.
### Moment of Inertia of a Disk Formula
  So the moment of inertia for a homogeneous disk about its center becomes

$$
    I = \frac{M}{\pi R^2} \int_0^{2 \pi} \int_0^R r^3 dr d\theta = \frac{2M}{R^2} \int_0^R r^3 dr
$$
$$
 = \frac{2M}{R^2} \Bigg[ \frac{r^4}{4} \Bigg]_0^R = \frac{2M}{R^2} \frac{R^4}{4} = M \frac{R^2}{2}.
$$

### Parallel Axis Theorem on Disk
 Now, what happens if the axis of rotation was shifted (but still parallel to the axis through the center of mass)? Well luckily for us we don't have to do a super complicated integral, we can just use the parallel axis thereom which says

$$
    I = I_\text{CM} + Md^2,
$$

where $I_\text{CM}$ is the moment of inertia through the center of mass, in this case $M \frac{R^2}{2}$. $d$ is the distance between the axis of rotation and the rotational axis throught the center of mass. You can test this formula by calculating the moment of inertia at an end of the rod. In that case we have $I_\text{CM} = M \frac{L^2}{12}$ and $d=\frac{L}{2}$.

 ### Moment of Inertia of a Disk With a Hole
 So what if the disk/cylinder is hollow? Then we simply change the limits of integration and the density. If the disk is hollow up to $r=a$ and then has mass up to $r=b$ the 2D density changes to

$$
    \frac{M}{\pi (b^2-a^2)}.
$$

So we get the following moment of inertia

$$
    I = \frac{2M}{b^2-a^2} \int_a^b r^3 dr = \frac{2M}{b^2-a^2} \Bigg[ \frac{r^4}{4} \Bigg]_a^b
$$
$$
 = \frac{2M}{b^2-a^2} \frac{b^4-a^4}{4} = M \frac{b^2+a^2}{2}.
$$

### Moment of Inertia of Cylindrical Shell
 Here we can see that if the disk/cylinder only consisted of a shell ($a \rightarrow b$) we get

$$
    I = M b^2.
$$

So now you know how to calculate the moment of inertia of a rod and a disk.

The Moment of Inertia of a Rod and a Disk

Does the Golden Gate Bridge's suspension wire curve like a parabola or catenary? Let's explore the mechanics behind these cable shapes.

**What makes the suspension wire of the Golden Gate Bridge curve like a parabola? Let us get into the mechanics of the shape of the suspension bridge's cables and what determines them to curve like a parabola or a catenary.**

## The Assumptions about the Problem

Some assumptions about the problem have to be made first. The main wire that the bridge rests on is assumed to be inextensible, homogeneous and infinitely flexible. The bridge is also assumed to be homogeneous. The main wire is connected at each end of the bridge. The wires that ties together the bridge and the main wire is assumed to be in an infinite amount so the bridge's load function is constant. The origin of coordinate system is located at the bottom of the main wire.

## Equilibrium of an Infinitesimal Wire Element

Lets make an free body diagram of an infinitesmal wire element of the main wire:

![Image of free body diagram of an infinitesmal wire element](https://i.imgur.com/vGGtqIv.png "Free body diagram of an infinitesmal wire element")

Here $\vec{T}$ is the tension of the cable, $d\vec{T}$ the infinitesmal change in the tension, $ds$ the infinitesmal length of the wire element and $dx$ the length in the $x$-direction. $q_c$ and $q_b$ is respectively the cable's and bridge's linear mass density. The equilibrium of the forces is

$$
    \vec{T} + d\vec{T} - \vec{T} - (q_c ds + q_b dx)g \hat{y} = 0.
$$

With some cancellation we get

$$
    d\vec{T} = (q_c ds + q_b dx)g \hat{y}.
$$

Lastly, the Pythagorean theorem gives us 

$$
    ds = \sqrt{(dx)^2 + (dy)^2} = \sqrt{1 + (y')^2} dx.
$$

So we get

$$
    d\vec{T} = (q_c \sqrt{1 + (y')^2} dx + q_b dx)g \hat{y}.
$$

## The ODE which Dictates the Shape of The Wire

The tension at the bottom of the main wire has only a $x$-component which we'll call $T_0$. So $\vec{T}(0) = T_0 \hat{x}$. If we do some integration of the equilibrium equation from earlier we get

$$
    \vec{T} = T_0 \hat{x} + \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right)g \hat{y},
$$

here we are just using $\tilde{x}$ in the integral so the upper limit, $x$, of the integral is correct. Also, as you might have noticed, the integral is the equation for the total length of a curve $y$ (<a href="https://www.amazon.com/Cable-Supported-Bridges-Concept-Design/dp/0470666285?crid=MH6235C2889K&dib=eyJ2IjoiMSJ9.dxWbL7XqIrrTulEQxoCUMmXNtIBKBNOGSJhIcZOZGBfweTsdihigyHeBfPV66CmOarA06iJfAuDhV23snJ7voMekG15MwV0o_iTCs5YdBAsJdy6QE6-frfED7i4Tlkb274YhoVZ49Jf4ek-aerMDJbXQ2Tj1QeKh22uPSpF8Z1lrwm1Ov5x0OX4y-lTV9gv8w5cM1R89HJvVe7Ck8YCuDo8lK_Q4QdIDGspukPsl8tI.NqSqY00VTQkh1G6T7FbKG1FSTevfln8tIW39dtSdFR0&dib_tag=se&keywords=cable+supported+bridges&qid=1734113319&sprefix=cable+support%2Caps%2C251&sr=8-1&linkCode=ll1&tag=engui-20&linkId=45c7ffd93fd045d6efbac260142fc9a6&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 133</a>). In the $x$-direction we now have

$$
    T \cos{\theta} = T \frac{dx}{\sqrt{(dx)^2 + (dy)^2}} = T_0,
$$

where $\theta$ is the angle between $\vec{T}$ and the $x$-axis. So

$$
    T = T_0 \frac{\sqrt{(dx)^2 + (dy)^2}}{dx}.
$$

In the $y$-direction we have

$$
    T \sin{\theta} = T \frac{dy}{\sqrt{(dx)^2 + (dy)^2}}
$$
$$
= T_0 \frac{dy}{dx} = \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right)g.
$$

We now have an expression for $y'$:

$$
    y' = \left(q_c \int_0^x \sqrt{1 + y'(\tilde{x})^2} d\tilde{x} + q_b x \right) \frac{g}{T_0}
$$

Lets do an additional derivation to get rid of the integral:

$$
    y'' = \frac{q_b g}{T_0} + \frac{q_c g \sqrt{1+(y')^2}}{T_0}.
$$

This is the differential equation that dictates the curve of suspension bridge.

## Is The Golden Gate Bridge Shaped as a Parabola or Catenary?

So what shape is the Golden Gate Bridge? To answer this lets look at two extremes, neglecting either the cable's mass or the bridge's mass. If we neglect the cable's load function i.e. its mass we get the differential equation

$$
    y'' = \frac{q_b g}{T_0}.
$$

With some simple integration we get

$$
    \frac{q_b g}{2 T_0} x^2 + bx + c,
$$

where $b$ and $c$ are constants. The origin is loacted at the bottom of the cable so $c=0$ and the cable is symmetric so $b=0$. Here we have our parabola. Neglecting the cable's mass is a pretty reasonable approximation since the bridge's mass is likely much larger than the cable's mass. Can you guess what we get if we do the opposite? If we neglect the bridge's mass we get the following differential equation:

$$
    y'' = \frac{q_c g \sqrt{1+(y')^2}}{T_0}.
$$

Since we already have pretty good hunch that the solution is in the shape of a catenary, lets check if $y = a \cosh{\frac{x}{a}}$ works. Then the left side is

$$
    y'' = \frac{\cosh{\frac{x}{a}}}{a}
$$

and the right side is

$$
    \frac{q_c g \sqrt{1+(y')^2}}{T_0} = \frac{q_c g \sqrt{1+\sinh^2{\frac{x}{a}}}}{T_0} =  \frac{q_c g \cosh{\frac{x}{a}}}{T_0}.
$$

For this to work $a$ needs to be

$$
    a = \frac{T_0}{q_c g}.
$$

That is the solution if the bridge's mass can be neglected. I don't know what sort of bridge that is... But hey now you know. I'll leave the harder task of solving the ODE with both masses present for you.

Is the Suspension Bridge Shaped like a Parabola or a Catenary?

To study transient processes in circuits, we use the Laplace transform. Here's how to transform a resistor, capacitor, and inductor to the s-domain.

**When studying transient processes, like when a switch is changed, in electrical circuits we have to make use of the Laplace transform. Here is a quick rundown on how to Laplace transform the resistor, capacitor and inductor to the s-domain.**

## Laplace Transform of Voltage and Current

First we need to establish what happens when we Laplace transform the current and voltage. We can see the voltage and current as "unknown" functions dependent on time, therefore their Laplace transforms become:

$$
        \begin{aligned}
            & i(t) \Rightarrow I(s) \\
            & u(t) \Rightarrow U(s)
        \end{aligned}
        
$$

 (<a href="https://www.amazon.com/Circuit-Analysis-Dummies-John-Santiago/dp/1118493125?crid=217JYBLWPN4MU&dib=eyJ2IjoiMSJ9.u-1JiLJprxLT2c6jXVAzjJSrTFPwJ4SIUBRMNUJWvfUbLeOVJokbDJHalMrmmDLFKqf9zuC7zaicNgmURouLIT5ZBsnmoX6VlHNa4gKoDo6J0g4rsBLiiCVzuyK2KJxjde-2jYVS4WhCwQMTo_LDoe3_turh8by37qBoEUZVpnxx5WcvJ9s-x1y_fJvjf-H4f7AXu2x4eX2lUfyT4sA31QqMYcDBkuaug3oDPsjcM2s.jniDHePI7W2mihyj3HYxKux8yRpBbTWKOQoNdkzAY98&dib_tag=se&keywords=circuit+analysis+for+dummies&qid=1733582045&sprefix=circuit+analysis%2Caps%2C192&sr=8-1&linkCode=ll1&tag=engui-20&linkId=543ba92ec6ba5c987774228967d40d8e&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 297</a>).

## Laplace Transform of the Resistor

So lets start with the easy component: the resistor. We know from Ohm's law that:

$$
    u(t) = R i(t).
$$

If we Laplace transform this we get:

$$
    U(s) = R I(s),
$$

since the resistor is time independent there is not much happening when moving from the t-domain to the s-domain (<a href="https://www.amazon.com/Circuit-Analysis-Dummies-John-Santiago/dp/1118493125?crid=217JYBLWPN4MU&dib=eyJ2IjoiMSJ9.u-1JiLJprxLT2c6jXVAzjJSrTFPwJ4SIUBRMNUJWvfUbLeOVJokbDJHalMrmmDLFKqf9zuC7zaicNgmURouLIT5ZBsnmoX6VlHNa4gKoDo6J0g4rsBLiiCVzuyK2KJxjde-2jYVS4WhCwQMTo_LDoe3_turh8by37qBoEUZVpnxx5WcvJ9s-x1y_fJvjf-H4f7AXu2x4eX2lUfyT4sA31QqMYcDBkuaug3oDPsjcM2s.jniDHePI7W2mihyj3HYxKux8yRpBbTWKOQoNdkzAY98&dib_tag=se&keywords=circuit+analysis+for+dummies&qid=1733582045&sprefix=circuit+analysis%2Caps%2C192&sr=8-1&linkCode=ll1&tag=engui-20&linkId=543ba92ec6ba5c987774228967d40d8e&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 297</a>). The resistor remains the same in both domains.

## Laplace Transform of a Derivative

The inductor and capacitor is time dependent, so they get a little bit more trickier. First we need to see what happens when we apply the Laplace transform to a derivative of a function $f(t)$. The Laplace transform of an arbitrary function $f(t)$ is given by:

$$
    F(s) = \int_{0-}^\infty e^{-st} f(t) dt
$$

(<a href="https://www.amazon.com/Circuit-Analysis-Dummies-John-Santiago/dp/1118493125?crid=217JYBLWPN4MU&dib=eyJ2IjoiMSJ9.u-1JiLJprxLT2c6jXVAzjJSrTFPwJ4SIUBRMNUJWvfUbLeOVJokbDJHalMrmmDLFKqf9zuC7zaicNgmURouLIT5ZBsnmoX6VlHNa4gKoDo6J0g4rsBLiiCVzuyK2KJxjde-2jYVS4WhCwQMTo_LDoe3_turh8by37qBoEUZVpnxx5WcvJ9s-x1y_fJvjf-H4f7AXu2x4eX2lUfyT4sA31QqMYcDBkuaug3oDPsjcM2s.jniDHePI7W2mihyj3HYxKux8yRpBbTWKOQoNdkzAY98&dib_tag=se&keywords=circuit+analysis+for+dummies&qid=1733582045&sprefix=circuit+analysis%2Caps%2C192&sr=8-1&linkCode=ll1&tag=engui-20&linkId=543ba92ec6ba5c987774228967d40d8e&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 274</a>). The $0-$ can be looked at more or less as what happened exactly before significant change happened to the function $f(t)$ (we'll get back to this). So lets preform the Laplace transform on $f'(t)$:

$$
    \int_{0-}^\infty e^{-st} f'(t) dt.
$$

With a little partial integration we get:

$$
    \bigg[e^{-st} f(t) \bigg]_{0-}^\infty - \int_{0-}^\infty -se^{-st} f(t) dt
$$
$$
= - f(0-) + sF(s) 
$$

(<a href="https://www.amazon.com/Circuit-Analysis-Dummies-John-Santiago/dp/1118493125?crid=217JYBLWPN4MU&dib=eyJ2IjoiMSJ9.u-1JiLJprxLT2c6jXVAzjJSrTFPwJ4SIUBRMNUJWvfUbLeOVJokbDJHalMrmmDLFKqf9zuC7zaicNgmURouLIT5ZBsnmoX6VlHNa4gKoDo6J0g4rsBLiiCVzuyK2KJxjde-2jYVS4WhCwQMTo_LDoe3_turh8by37qBoEUZVpnxx5WcvJ9s-x1y_fJvjf-H4f7AXu2x4eX2lUfyT4sA31QqMYcDBkuaug3oDPsjcM2s.jniDHePI7W2mihyj3HYxKux8yRpBbTWKOQoNdkzAY98&dib_tag=se&keywords=circuit+analysis+for+dummies&qid=1733582045&sprefix=circuit+analysis%2Caps%2C192&sr=8-1&linkCode=ll1&tag=engui-20&linkId=543ba92ec6ba5c987774228967d40d8e&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 276</a>).

 ## Laplace transform of the Inductor

With that out of the way lets preform the Laplace transform on the inductor first. We know that the potential over an inductor is given by:

$$
    u(t) = L \frac{d i(t)}{dt}.
$$

So with our knowledge of Laplace transforming a derivative we get the following:

$$
    U(s) = L (s I(s) - i(0-)).
$$

The inductor is said to be energized if there is an initial current $i(0-)$. If there is no initial current we get a simpler relation of:

$$
    U(s) = sL I(s).
$$

We can also rewrite the relation to get the current:

$$
    I(s) = \frac{U(s)}{sL} + \frac{i(0-)}{s}
$$

 (<a href="https://www.amazon.com/Circuit-Analysis-Dummies-John-Santiago/dp/1118493125?crid=217JYBLWPN4MU&dib=eyJ2IjoiMSJ9.u-1JiLJprxLT2c6jXVAzjJSrTFPwJ4SIUBRMNUJWvfUbLeOVJokbDJHalMrmmDLFKqf9zuC7zaicNgmURouLIT5ZBsnmoX6VlHNa4gKoDo6J0g4rsBLiiCVzuyK2KJxjde-2jYVS4WhCwQMTo_LDoe3_turh8by37qBoEUZVpnxx5WcvJ9s-x1y_fJvjf-H4f7AXu2x4eX2lUfyT4sA31QqMYcDBkuaug3oDPsjcM2s.jniDHePI7W2mihyj3HYxKux8yRpBbTWKOQoNdkzAY98&dib_tag=se&keywords=circuit+analysis+for+dummies&qid=1733582045&sprefix=circuit+analysis%2Caps%2C192&sr=8-1&linkCode=ll1&tag=engui-20&linkId=543ba92ec6ba5c987774228967d40d8e&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 298</a>). We can make circuit diagrams with these relations:

![Image of Laplace transform of the inductor](https://i.imgur.com/IkgvQUe.png "Laplace transform of the inductor")

Note the Thevenin/Norton equivalents under "Initial current".

## Laplace transform of the Capacitor

Lets move on to the capacitor. The procedure of Laplace transforming the capacitor is the same as the inductor. The current through a capacitor is given by:

$$
    i(t) = C \frac{du(t)}{dt}.
$$

So the Laplace transform of the current is:

$$
    I(s) = C(sU(s) - u(0-)).
$$

We can see here that if there is no initial voltage over the capacitor we get the simpler relation:

$$
    I(s) = sC U(s).
$$

Lastly the the Laplace transformed voltage is:

$$
    U(s) = \frac{I(s)}{sC} + \frac{u(0-)}{s}
$$

 (<a href="https://www.amazon.com/Circuit-Analysis-Dummies-John-Santiago/dp/1118493125?crid=217JYBLWPN4MU&dib=eyJ2IjoiMSJ9.u-1JiLJprxLT2c6jXVAzjJSrTFPwJ4SIUBRMNUJWvfUbLeOVJokbDJHalMrmmDLFKqf9zuC7zaicNgmURouLIT5ZBsnmoX6VlHNa4gKoDo6J0g4rsBLiiCVzuyK2KJxjde-2jYVS4WhCwQMTo_LDoe3_turh8by37qBoEUZVpnxx5WcvJ9s-x1y_fJvjf-H4f7AXu2x4eX2lUfyT4sA31QqMYcDBkuaug3oDPsjcM2s.jniDHePI7W2mihyj3HYxKux8yRpBbTWKOQoNdkzAY98&dib_tag=se&keywords=circuit+analysis+for+dummies&qid=1733582045&sprefix=circuit+analysis%2Caps%2C192&sr=8-1&linkCode=ll1&tag=engui-20&linkId=543ba92ec6ba5c987774228967d40d8e&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 298</a>). The circuit diagrams then become:

![Image of Laplace transform of the capacitor](https://i.imgur.com/za9EeTS.png "Laplace transform of the capacitor")

And that is it! You now know how to Laplace transform the resistor, inductor and capacitor!

Laplace Transform of DC Circuits to s-domain

The Navier-Stokes equation is unsolvable with turbulence, but solutions exist for laminar flow, like Couette flow or flow between two plates.

**The Navier Stokes equation is famously known for being unsolvable if there is any turbulence in a fluid but there are certain scenarios where laminar flow exist and a solution for the velocity profile can be derived. One of these scenarios is flow past two plates either caused by one plate moving, Couette flow, or a pressure gradient.**

Before solving these two problems we need to make some assumptions. First we assume that our fluid is incompressible and Newtonian so we can use the Navier Stokes equation and the mass continuity equation. Secondly, the flow has to be stationary (so all time dependent terms in the Navier Stokes equation is zero) and in one direction, the $x$-direction, so the velocity vector only has one component ($v=w=0$). The plates also has to be considered infinite. Lastly, the flow has to be laminar which has already been mentioned and the impact from gravity has to be considered negligible ($g = 0$).

## Couette Flow - One Plate Moves

So lets start with the case where one of the plates is moving. From our assumptions that fluid is Newtonian, stationary and goes in the $x$-direction we get the following mass continuity equation:

$$
    \frac{\partial u}{\partial x} + \cancel{\frac{\partial v}{\partial y}} + \cancel{\frac{\partial w}{\partial z}} = \frac{\partial u}{\partial x} = 0,
$$

where $u$, $v$ and $w$ is the velocity in the $x$-, $y$- and $z$-direction. It is then safe to assume that the velocity profile $u$ is only dependent on $y$ because the plates' areas are infinite. Now lets look at the Navier Stoke equation in all directions and see which terms we can cancel out because of our assumptions and the mass contuinity equation:

$$
    \left\{
        \begin{aligned}
            & \rho \left( \cancel{\frac{\partial u}{\partial t}} + \cancel{u \frac{\partial u}{\partial x}} + \cancel{v \frac{\partial u}{\partial y}} + \cancel{w \frac{\partial u}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial x}} + \mu \left(\cancel{\frac{\partial^2 u}{\partial x^2}} + \frac{\partial^2 u}{\partial y^2} + \cancel{\frac{\partial^2 u}{\partial z^2}} \right) + \cancel{\rho g_x} \\
            & \rho \left( \cancel{\frac{\partial v}{\partial t}} + \cancel{u \frac{\partial v}{ \partial x}} + \cancel{v \frac{\partial v}{\partial y}} + \cancel{w \frac{\partial v}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial y}} + \mu \left(\cancel{\frac{\partial^2 v}{\partial x^2}} + \cancel{\frac{\partial^2 v}{\partial y^2}} + \cancel{\frac{\partial^2 v}{\partial z^2}} \right) + \cancel{\rho g_y} \\
            & \rho \left( \cancel{\frac{\partial w}{\partial t}} + \cancel{u \frac{\partial w}{ \partial x}} + \cancel{v \frac{\partial w}{\partial y}} + \cancel{w \frac{\partial w}{\partial z}} \right) = \cancel{- \frac{\partial p}{\partial z}} + \mu \left(\cancel{\frac{\partial^2 w}{\partial x^2}} + \cancel{\frac{\partial^2 w}{\partial y^2}} + \cancel{\frac{\partial^2 w}{\partial z^2}} \right) + \cancel{\rho g_z}
        \end{aligned}
    \right. .
$$

As you can see, a bunch of terms cancel out and all that remains is just one single expression:

$$
    \mu \frac{\partial^2 u}{\partial y^2} = 0 \Rightarrow \frac{\partial^2 u}{\partial y^2} = 0.
$$

So if we integrate this expression twice we get:

$$
    u(y) = C_1 y + C_2,
$$

where $C_1$ and $C_2$ are constants. To determine these constants we have to use the no-slip boundary conditions of the plates. If we assume that the origin is located exactly in between the two plates we know that:

$$
    \left\{
        \begin{aligned}
            & u(+h) = u_p \\
            & u(-h) = 0
        \end{aligned}
    \right. ,
$$

where $h$ is the shortest distance from the origin to a plate and $u_p$ is the speed of the moving plate. This gives the constants:

$$
    \left\{
        \begin{aligned}
            & C_1 = \frac{u_p}{2 h} \\
            & C_2 = \frac{u_p}{2}
        \end{aligned}
    \right. .
$$

So our velocity profile becomes:

$$
    u(y) = \frac{u_p}{2} \left(1 + \frac{y}{h} \right).
$$

 (<a href="https://www.amazon.com/Fluid-Mechanics-Joseph-H-Spurk-ebook/dp/B0828MVGKG?crid=1D44S7EXN88JZ&dib=eyJ2IjoiMSJ9.SQ94Q1-GKZ7QuFcUVU-xi2UZoW_Fc1omG_OcDMaVNpZvaTSTS1Hriu_znh3-b5N-Z3D-DLd5pvVaeJ76I2Gs0ao1HAzbVbDZdzGLtwuILeGAzqwLDD5_QWCkDbCAAzpwx5A8TeXYtPnr7qIKpL2EgbKuxXoBHB_4w1-sBw0FLtQJr49XjF2inhBRYhyDQBpKVB3E0QjZBr3CND0tN-GeA_uHsLNveK6gNIZw_0WTGZQ.P0_SrAEDzPbrIAo_0CGCwFRJS3PRDIAcw9VF84dGJBA&dib_tag=se&keywords=fluid+mechanics+spurk&qid=1734115957&sprefix=fluid+mechanics+spurk%2Caps%2C265&sr=8-1&linkCode=ll1&tag=engui-20&linkId=98ce0c3d01a91ec433e60713dca30c91&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 180-181</a>). This solution is called the Couette flow which is named after the french physicist Maurice Couette. 

## Poiseuille Flow - Laminar Flow Past Two Plates Caused by Pressure Gradient

So what solution do we get if the two plates are stationary and we instead have a constant pressure gradient? The assumptions mentioned earlier are still the same so we get the same mass continuity equation. So lets see what cancels out in our Navier Stokes equations:

$$
    \left\{
        \begin{aligned}
            & \rho \left( \cancel{\frac{\partial u}{\partial t}} + \cancel{u \frac{\partial u}{\partial x}} + \cancel{v \frac{\partial u}{\partial y}} + \cancel{w \frac{\partial u}{\partial z}} \right) = - \frac{\partial p}{\partial x} + \mu \left(\cancel{\frac{\partial^2 u}{\partial x^2}} + \frac{\partial^2 u}{\partial y^2} + \cancel{\frac{\partial^2 u}{\partial z^2}} \right) + \cancel{\rho g_x} \\
            & \rho \left( \cancel{\frac{\partial v}{\partial t}} + \cancel{u \frac{\partial v}{ \partial x}} + \cancel{v \frac{\partial v}{\partial y}} + \cancel{w \frac{\partial v}{\partial z}} \right) = - \frac{\partial p}{\partial y} + \mu \left(\cancel{\frac{\partial^2 v}{\partial x^2}} + \cancel{\frac{\partial^2 v}{\partial y^2}} + \cancel{\frac{\partial^2 v}{\partial z^2}} \right) + \cancel{\rho g_y} \\
            & \rho \left( \cancel{\frac{\partial w}{\partial t}} + \cancel{u \frac{\partial w}{ \partial x}} + \cancel{v \frac{\partial w}{\partial y}} + \cancel{w \frac{\partial w}{\partial z}} \right) = - \frac{\partial p}{\partial z} + \mu \left(\cancel{\frac{\partial^2 w}{\partial x^2}} + \cancel{\frac{\partial^2 w}{\partial y^2}} + \cancel{\frac{\partial^2 w}{\partial z^2}} \right) + \cancel{\rho g_z}
        \end{aligned}
    \right. .
$$

We can see that the pressure gradient is zero in the $y$- and $z$-direction. So we get the following differential equation:

$$
    \mu \frac{\partial^2 u}{\partial y^2} = \frac{\partial p}{\partial x}.
$$

With some integration we get:

$$
    u(y) = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{y^2}{2} + C_1 y + C_2,
$$

where $C_1$ and $C_2$ are constants. The velocity at both plates are now zero because of the no-slip condition, so:

$$
    \left\{
        \begin{aligned}
            & u(+h) = 0 \\
            & u(-h) = 0
        \end{aligned}
    \right. .
$$

This gives the constants:

$$
    \left\{
        \begin{aligned}
            & C_1 = 0 \\
            & C_2 = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{h^2}{2}
        \end{aligned}
    \right. .
$$

So our velocity profile becomes:

$$
    u(y) = \frac{1}{\mu} \frac{\partial p}{\partial x} \frac{h^2}{2} \left(1 + \left(\frac{y}{h} \right)^2 \right),
$$

(<a href="https://www.amazon.com/Fluid-Mechanics-Joseph-H-Spurk-ebook/dp/B0828MVGKG?crid=1D44S7EXN88JZ&dib=eyJ2IjoiMSJ9.SQ94Q1-GKZ7QuFcUVU-xi2UZoW_Fc1omG_OcDMaVNpZvaTSTS1Hriu_znh3-b5N-Z3D-DLd5pvVaeJ76I2Gs0ao1HAzbVbDZdzGLtwuILeGAzqwLDD5_QWCkDbCAAzpwx5A8TeXYtPnr7qIKpL2EgbKuxXoBHB_4w1-sBw0FLtQJr49XjF2inhBRYhyDQBpKVB3E0QjZBr3CND0tN-GeA_uHsLNveK6gNIZw_0WTGZQ.P0_SrAEDzPbrIAo_0CGCwFRJS3PRDIAcw9VF84dGJBA&dib_tag=se&keywords=fluid+mechanics+spurk&qid=1734115957&sprefix=fluid+mechanics+spurk%2Caps%2C265&sr=8-1&linkCode=ll1&tag=engui-20&linkId=98ce0c3d01a91ec433e60713dca30c91&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">p. 181-183</a>). The velocity profile is a parabola with its maximum at the middle of the plates. 

![Pressure Gradient Flow Between Two Plates](https://i.imgur.com/xIiH6ge.png "Pressure Gradient Flow Between Two Plates")

So that is it, that is how you calculate fluid flow past two plates!

Couette and Poiseuille Flow - Laminar Flow Past Two Plates

The chemical potential varies with temperature. Let's explore its calculation and approximations in an n-type semiconductor.

**The chemical potential has different characteristics depending on the temperature. Lets calculate and understand the approximations of the chemical potential in an electron doped semiconductor.**

## The Chemical Potential at Low Temperatures

Lets start with temperatures close to absolute zero. You have probably heard or read that the chemical potential must lie in the center between the donor level and the bottom of the conduction band because of charge neutrality. Lets do some math to find out why that is.

### Bandgap Energy, Donor Energy, and Donor Level

Lets use the variable $E_g$ for the bandgap energy, $E_d$ is the donor energy so the the donor level is situated at $E_g-E_d$. $\mu$ is our desired chemical potential. $n$ is the electron density in the conduction band, $n_d$ is the donor density and $n_d^+$ is the density of the ionized donors. From the Fermi-Dirac distribution we know that the probability of an electron at the energy $E_g-E_d$ and the temperature $T$ is:

$$
    f(E_g - E_d) = \frac{1}{e^{(E_g-E_d-\mu)/(k_B T)}+1}.
$$

So the donor electron density at the donor level is:

$$
    \frac{n_d}{e^{(E_g-E_d-\mu)/(k_B T)}+1}.
$$

Every donor electron not at the donor level is at the conduction band so:

$$
    n_d^+ = n_d (1 - \frac{1}{e^{(E_g-E_d-\mu)/(k_B T)}+1}) 
$$
$$
= \frac{n_d}{e^{-(E_g-E_d-\mu)/(k_B T)}+1}.
$$

Since we are at low temperatures $T \ll E_d/k_B$. This means $-(E_g-E_d-\mu)/(k_B T) \gg 1$ so we can make the following approximation (the 1 in the denominator is neglected):

$$
    n_d^+ \approx n_d e^{(E_g-E_d-\mu)/(k_B T)}.
$$

Having low temperatures also means that the electrons at the valence band don\'t have enough energy to make an excitation to the conduction band meaning $n \approx n_d^+$. But the electron density $n$ can also be expressed as:

$$
    n = \frac{1}{\sqrt{2}} (\frac{m_e^* k_B T}{\pi \hbar^2})^{3/2} e^{-(E_g-\mu)/(k_B T)} 
$$
$$
= n_0 e^{-(E_g-\mu)/(k_B T)}
$$

here is $m_e^*$ the effective mass of the electrons in the conduction band. Lets put the different expressions for $n$ at low temperature equal to each other to get the chemical potential:

$$
    n_d e^{(E_g-E_d-\mu)/(k_B T)} = n_0 e^{-(E_g-\mu)/(k_B T)}
$$

So now we get the chemical potential:

$$
    \mu = E_g - \frac{1}{2} E_d - \frac{1}{2} k_B T \ln{(\frac{n_0}{n_d})}.
$$

So when the temperature approaches zero, $T \to 0$, the chemical potential becomes:

$$
    \mu = E_g - \frac{1}{2}E_d,
$$

since polynomial functions dominates over logarithmic functions ($n_0$ contains $T^{3/2}$) when it comes to limits.

## The Chemical Potential at Medium Temperatures

So that is how the chemical potential behaves at low temperatures, how about "medium" temperatures like room temperature? Well here we have more thermal energy to work with so pretty much every donor atom is ionized ($n_d^+ \approx n_d$) but the thermal energy is still not big enough to excite a significant amount of electrons in the valence band so $n \approx n_d$. This makes the math easier, we now have:

$$
    n = n_0e^{-(E_g-\mu)/(k_B T)} = n_d,
$$

so:

$$
    \mu = E_g + k_B T \ln{(\frac{n_d}{n_0})},
$$

for "medium" temperatures.

## The Chemical Potential at High Temperatures

For high temperatures, like 1200 K, electron holes has manifested at the valence band and now plays the bigger role in the electron density in the conduction band. So now $n \approx p$, where $p$ is the electron hole density at the valence band.

### Intrinsic State

 The semiconductor is said to be in an intrinsic state, meaning that it can be seen as undoped. So now we have the following relation:

$$
    \frac{1}{\sqrt{2}} (\frac{m_e^* k_B T}{\pi \hbar^2})^{3/2} e^{-(E_g-\mu)/(k_B T)} 
$$
$$
= \frac{1}{\sqrt{2}} (\frac{m_h^* k_B T}{\pi \hbar^2})^{3/2} e^{-\mu/(k_B T)},
$$

where $m_h^*$ is the effective mass of the electron holes. This leads to the chemical potential of:

$$
    \mu = \frac{1}{2} E_g + \frac{3}{4} k_B T \ln{(\frac{m_h^*}{m_e^*})},
$$

(<a href="https://www.amazon.com/Solid-State-Physics-Philip-Hofmann/dp/352741410X?crid=9TFN7QN1X3JZ&dib=eyJ2IjoiMSJ9.tThGMryvCmAutjgHm3s1vIreko5fGPWxa-iRNtwqxrrDbBefHOjqdOvkpTPh2OJjjfWEZ7c522o5JrziznhHbQMj6Axs4LsPdH-gJKzK1jSdOoagscymyNfTSYGhk7HNTjr3_qEspf07WaTzNm95dd0HazsKH3p1czaHMxPFgeGLOUXUltEA-haeuoVgiPSBzoK3Vy0b97kyCyzDT1wc-M3mTkvPIzuWL_W_0Yx05tA.-Sv304NOx4dyCzcKczRUqai6ckbfm1W4FFYTUZsdbK0&dib_tag=se&keywords=solid+state+physics+an+introduction&qid=1734117174&sprefix=solid+state+physics+an+introduction%2Caps%2C218&sr=8-3&linkCode=ll1&tag=engui-20&linkId=a05b80ea8a819ceedc7a3fc3d6d23302&language=en_US&ref_=as_li_ss_tl" rel="nofollow sponsored">eq 7.15</a>). The chemical potential is located midway between the conduction band and the valence band if $m_e^* = m_h^*$.

So there we have it! We have now calculated the chemical potential for different temperature regions for an n-type semiconductor, although we did it with some significant approximations but without those this would get a lot harder if not impossible.

The Chemical Potential of an n-type Semiconductor

Why is it that a falling chimney breaks at a third of its length? Lets solve this physics problem!

**The falling chimney problem is the perfect mix of mechanical dynamics and beam theory. By just using the fundamentals of those two fields we can accurately predict where the chimney cracks when falling, under some assumptions...**

So lets get the assumptions of the falling chimney problem out of the way first. The first assumption is that there is no air resistance. If air resistance had been present the top of the chimney would experience a greater stress from the air than the lower part which makes our calculations much more complex and this is not an article about fluid mechanics. The second assumption is that the chimney is homogeneous, meaning that it has the same diameter across its length and the same material. If it had been inhomogeneous certain parts would handle stress better than others.

## Angle Acceleration of the Chimney

Now when that is out of the way, lets make some calculations! The first thing we'll have to determine is the angle acceleration $\ddot{\theta}$. To do that we use a beam's, in this case a chimney's, moment of inertia when attached at one of its ends:

$$
    I = \frac{1}{3} m L^2.
$$

Here $m$ is the chimneys mass and $L$ is its length. The moment from the gravitation acting on the chimney must equal that of the angle acceleration:

$$
    mg \frac{L}{2} \sin{\theta} = I \ddot{\theta} = \frac{1}{3}mL^2 \ddot{\theta}.
$$

So now we have our angle acceleration $\ddot{\theta}$ and apparently it is dependent on $\theta$:

$$
    \ddot{\theta} = \frac{3g}{2L} \sin{\theta}.
$$

Is this plausible? Well $\ddot{\theta}$ increases with $g$ and decreases with $L$ which makes sense. At $\theta=0$ the angle acceleration is zero which also makes sense since at that angle the chimney is "standing up".

## The Chimney's Load Function

If the chimney was locked and stationary at an angle $\theta$ the load function would be:

$$
    w(d) = -\frac{mg}{L} \sin{\theta}.
$$

If we go to the other extreme and say that the chimney was in free fall it would lead to the load function being zero. So falling negates the load function and in our case the chimney is rotating. So $\ddot{\theta}$ has a similar effect to free fall. So the load function in our case will be the following:

$$
    w(d)=\frac{dm}{L} \ddot{\theta} - \frac{mg}{L} \sin{\theta} 
$$
$$
= \frac{3dmg}{2L^2} \sin{\theta} - \frac{mg}{L} \sin{\theta}.
$$

We can see here that the load function is zero at $d = \frac{2L}{3}$. This makes sense because $\ddot{\theta}$ is equal to $g \sin{\theta}$ at $d = \frac{2L}{3}$. In other words, this point is in a "virtual" free fall.

## The Chimney's Shear Force Function

The relation between the load function and the shear force function is:

$$
    -w(d)=V'(d).
$$

If we integrate that we get:

$$
    V(d) = -\frac{3d^2mg}{4L^2} \sin{\theta}+ \frac{dmg}{L} \sin{\theta}+ C,
$$

where $C$ is a constant. We can determine this constant with the boundary condition $V(L) = 0$. In other words: the shear force is zero at the end of the chimney which is logical. This means:

$$
    C = -\frac{mg}{4}.
$$

If we take the absolute value of the shear force function we can check where the force is the greatest:

$$
    |V(d)|=\left|mg\, \sin{\theta} \left(-\frac{3d^2}{4L^2} + \frac{d}{L} - \frac{1}{4}\right) \right|.
$$

Here we can see that the shear force is the greatest when $d \to 0$ on the interval $0< d < L$, so if the shear force was the cause of the crack in the chimney it would be at the base. But that is not the only place where the chimney can break... 

## The Chimney's Bending Moment Function

The final step is looking at the bending moment function which has the following relation with the shear force function:

$$
    M'(d)=V(d).
$$

So with a bit of integration we get:

$$
    M(d)=mg\, \sin{\theta} \left(-\frac{d^3}{4L^2} + \frac{d^2}{2L} - \frac{d}{4}\right) + K,
$$

where $K$ is a constant. We'll once again use the boundary condition at the end of the chimney, which in this case is also zero. This leads to $K=0$. So where is $|M(d)|$ the greatest value? At $d=0$ the bending moment is zero so it is not there. Well lets check where its derivative, the shear force function, is zero to get the maximums:

$$
    M'(d)=V(d)=0 \Longrightarrow -\frac{3d^2}{4L^2} + \frac{d}{L} - \frac{1}{4} = 0.
$$

From this we get the roots $d=L$ and $d=\frac{L}{3}$. $d=L$ does not make much sense since its the end of the chimney, you cant real chop of something that isn't there... So we get $d=\frac{L}{3}$ and that is why a falling chimney breaks at a third of its length from its base.

<div style="position: relative;
    padding-bottom: 56.25%; /* 16:9 */
    height: 0;">
<iframe style="position: absolute;
    top: 0;
    left: 0;
    width: 100%;
    height: 100%;" width="560" height="315" src="https://www.youtube.com/embed/efjegZ6aWKI?start=149" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

</div>

The video above shows a demolition of a smokestack, as you can see it breaks close to a third of its length when it falls.

The Falling Chimney Physics Problem

Calculating the energy of a two-electron stationary state is complex, but methods like perturbation theory and the variational method offer good estimates.

**Calculating the energy of the stationary state of a two electron system, like a helium atom or lithium ion, gets quiet a bit more tricky than for a single electron system. In fact, calculating the Schrödinger equation is impossible. Luckily we have methods like perturbation theory and the variational method to estimate the energy.**

First we need to establish the Hamiltonian for the system which is:
$$
    \hat{H} = -\frac{\hbar}{2 m_e} \nabla^2_1 - \frac{\hbar}{2 m_e} \nabla^2_2 - \frac{Ze^2}{4 \pi \epsilon_0 r_1 } - \frac{Ze^2}{4 \pi \epsilon_0 r_2 } + \frac{e^2}{4 \pi \epsilon_0 | \vec{r}_1 - \vec{r}_2 | }.
$$
$\frac{\hbar}{2 m_e} \nabla^2_1$ is the kinetic energy for the first electron, same goes for the second. $\frac{Ze^2}{4 \pi \epsilon_0 r_1 }$ is Coulomb's law between the atom core's charge $Ze$ and the first electron's charge $e$. Lastly we have $\frac{e^2}{4 \pi \epsilon_0 | \vec{r}_1 - \vec{r}_2 | }$ which is Coulomb's law between the two electrons, this term is positive because the electrons repel each other.

## Stationary State Estimation Two Independent Electrons

The first estimation we'll do is having the two electrons not react with each other which means that the Coulomb's law term between them is zero. Then the combined wave function of the system is the following:

$$
    \psi_+ (\vec{r}_1, \vec{r}_2) = \psi_{100}(r_1) \psi_{100}(r_2),
$$

where

$$
    \psi_{100}(r) = \frac{Z^{3/2}}{\sqrt{\pi} a_0^{3/2}} e^{Zr/a_0}
$$

is the wave function of a single electron system, $a_0$ is the Bohr radius. So our wave function of the two electron system is

$$
    \psi_+ (r_1, r_2) = \frac{Z^{3}}{\sqrt{\pi} a_0^{3}} e^{Z(r_1+r_2)/a_0},
$$

which leads the stationary state energy of:

$$
    E = 2Z^2 \cdot (-13.6 \text{ eV}).
$$

## Stationary State Estimation Using Perturbation Theory

But that is not so satisfactory nor a good estimate, which leads to our second estimation method: perturbation theory. Lets use the first order correction to estimate the impact that Coulomb's law between the two electrons has on the energy:

$$
    E^{(1)} = \langle \psi_+ | \frac{e^2}{4 \pi \epsilon_0} \frac{1}{|\vec{r}_1 - \vec{r}_2|} | \psi_+ \rangle = \frac{5}{4} Z \cdot 13.6 \text{ eV}.
$$

So now our estimated stationary state energy is:

$$
    E = 2Z^2 \cdot (-13.6 \text{ eV}) + \frac{5}{4} Z \cdot 13.6 \text{ eV}.
$$

## Stationary State Estimation Using Variational Method

Lets try getting an even better estimation with the variational method. If we change the atom number $Z$ to a parameter $\alpha$ we get:

$$
    \psi_0 (r_1, r_2) = \frac{\alpha^{3}}{\sqrt{\pi} a_0^{3}} e^{\alpha (r_1+r_2)/a_0}.
$$

We now use this wave function as the base to calculate the expectation value of the Hamiltonian:

$$
    E(\alpha) = \langle \psi_0 (r_1, r_2) | \hat{H} | \psi_0 (r_1, r_2) \rangle 
$$
$$
 = (-2 \alpha^2 + 4 Z \alpha - \frac{5}{4} \alpha) \cdot (-13.6 \text{ eV}).
$$

Here we can see that if we choose $\alpha = Z$ we get the same result as for the first order correction perturbation theory. The variational method says $E_\text{ground} \leq E(\alpha)$ so lets check where the minimum of $E(\alpha)$ is:

$$
    \frac{\partial E(\alpha)}{\partial \alpha} = (-4 \alpha + 4Z - \frac{5}{4}) = 0 \Longrightarrow \alpha = Z - \frac{5}{16}.
$$

So our estimated stationary state energy with the variational method is

$$
    E_\text{ground} = (Z-\frac{5}{16})^2 \cdot (-13.6 \text{ eV}).
$$

## Helium Atom Variational Method

For example, the estimated stationary state energy with the variational method for helium ($Z=2$) is $-77.5$ eV which is close to the experimental value of $-78.975$ eV, in other words this method of estimation works very well! You now know a couple of methods to estimate the stationary state energy for electrons!

FILTERED POSTS

Benford's Law - Analysis and Example

Momentum Equation - Linear And Angular Conservation

Continuity Equation for Fluids in Different Coordinate Systems

Undamped and Damped Oscillations in Free Vibration

Relative Motion, Velocity, and Acceleration

RECENT POSTS