Dear users

I usually write my perturbations functions in python to see the result that i want and then i write them in the ASPECT muparser format.

For example, i wanted to edit this CITCOM function perturbation that i found in a paper:

subsection Function

set Coordinate system = spherical

set Variable names = r,theta

```
set Function constants = r_inner=3481000, r_outer=6371000, p=300.0, frequency=6, w=2
set Function expression = p*sin(pi*(r-r_inner)/(r_outer-r_inner))*sin(frequency*theta/w)
```

end

And then adding small and large perturbations. So i wrote this python code:

```
import numpy as np
import matplotlib.pyplot as plt
# Parameters
r_inner = 3481000 # Inner radius in meters
r_outer = 6371000 # Outer radius in meters
p_large = 700.0 # Amplitude of large perturbations
p_small = 100.0 # Amplitude of small perturbations
frequency_large = 3 # Frequency of large perturbations
frequency_small = 30 # Frequency of small perturbations
spread_angle_large_degrees = 90 # Spread angle of large perturbations in degrees
start_angle_degrees = 270 # Starting angle for the perturbation in degrees
# Convert angles from degrees to radians
spread_angle_large = np.radians(spread_angle_large_degrees)
start_angle = np.radians(start_angle_degrees)
# Temperature function
def temperature_function(r, theta):
# Adjust theta to shift the start angle of the perturbation
adjusted_theta = (theta + start_angle) % (2 * np.pi)
# Determine the large perturbation by using a step function approach
if (adjusted_theta < spread_angle_large) or (adjusted_theta > 2 * np.pi - spread_angle_large):
# Inside the large perturbation window
temp = p_large * np.sin(np.pi * (r - r_inner) / (r_outer - r_inner)) * np.sin(frequency_large * theta)
else:
# Outside the large perturbation window, apply small perturbations
temp = p_small * np.sin(np.pi * (r - r_inner) / (r_outer - r_inner)) * np.sin(frequency_small * theta)
return temp
# Grid for r and theta values
r_values = np.linspace(r_inner, r_outer, 100)
theta_values = np.linspace(0, 2 * np.pi, 100)
r_grid, theta_grid = np.meshgrid(r_values, theta_values)
# Calculate temperature values on the grid
temperature_grid = np.vectorize(temperature_function)(r_grid, theta_grid)
# Plotting
plt.figure(figsize=(12, 6))
X, Y = r_grid * np.cos(theta_grid), r_grid * np.sin(theta_grid)
plt.contourf(X, Y, temperature_grid, 100, cmap='viridis')
plt.colorbar()
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Adjusted Temperature Distribution in a Spherical Shell (2D) with Starting Angle')
plt.show()
```

To get this:

So i get 3 large perturbations and a lot of small perturbations. The problem is that I tried the whole day to put this expression into ASPECT with poor results, this is the best i could do because of the complexity:

#subsection Function

# set Coordinate system = spherical

#set Variable names = r,theta

```
set Function constants = r_inner=3481000, r_outer=6371000, p_large=700.0, p_small=100.0, frequency_large=3*2*pi, frequency_small=30*2*pi, spread_angle_large=180*pi/180, start_angle=270*pi/180, pi=3.14159265359
set Function expression = (0.5*(1+sign(sin(frequency_large*(theta+start_angle)/2)))*sign(sin(spread_angle_large*frequency_large/2-sin(frequency_large*(theta+start_angle)/2))))*p_large*sin(pi*(r-r_inner)/(r_outer-r_inner))*sin(frequency_large*theta) + (1-(0.5*(1+sign(sin(frequency_large*(theta+start_angle)/2)))*sign(sin(spread_angle_large*frequency_large/2-sin(frequency_large*(theta+start_angle)/2)))))*p_small*sin(pi*(r-r_inner)/(r_outer-r_inner))*sin(frequency_small*theta)
end
```

The function is very bugged and is not clean like my python code.

Can someone please help me to convert this function?

Thank you in advance for the help

Francesco