Efield and precipitation input generation handled incongruently

[mattzett]

What happened?

In the creation of electric field input data, the time dependence of the field are handled through the user-specified function; whereas for precipitation they are handled from the pygemini function which prevents users from easily giving both the shape and time-dependence of the precipitation they want without modifying pygemini itself.

Relevant log output

No response

[mattzett]

Upon further review, we appear to have an example here: https://github.com/gemini3d/gemini-examples/tree/main/init/aurora_EISCAT3D which has a user-defined precipitation that varies in time.

[mattzett]

We probably need some readmes to cover the most common use cases; I frequently forget the extent to which we've included these in the examples and I think there are so many examples that it can be difficult to find the most pertinent ones. We may want to consider separating them according to "curated" vs. "old but potentially useful"

[jdiazpena]

I did something similar to have precipitation varying in time and space for the SAID simulations, though it is hard coded. I edited the particles_BCs and Efiel_BCs to accept functions for both Qprecip and Eprecip (as well as the Efield, which was already possible). For python I think it is the init on each folder.

I have to come up with a gaussian or standard function for the Eprecip if a function is not given, but it was done for my simulations at least.

if isfield(p, "Qprecip_function")
Qfunc = str2func(p.Qprecip_function);
else
Qfunc = str2func("gemini3d.particles.gaussian2d");
end

if isfield(p, "E0precip_function")
Efunc = str2func(p.E0precip_function);
end

for i = i_on:i_off
precip.Qit(:,:,i) = Qfunc(precip,p.Qprecip,p.Qprecip_background);
% precip.E0it(:,:,i) = p.E0precip;
precip.E0it(:,:,i) = Efunc(precip,p.E0precip,p.E0precip_background);
end

[mattzett]

Yeah I was initially thinking it would have to be hard-coded in the the mat_gemini or pygemini scripts, but I'm not actually sure that is the case because of that one example above which appears to just have it in a user-defined area. I'm actually running that example again to see it it works...

[jdiazpena]

Yes, I looked at that example to build mine. Pygemini and mat_gemini are currently coded to take functions for Q, but not for E0... I added the E0 part to it, not hard coded in the sense that it is in pygemini what values I use, but that I hardcoded that you always have to give a function of E0 (I did not creat a default function of E0, while a default function of Q does exist, which is a gaussian)

[mattzett]

Ah I see, so the issue is really that we need to allow for user-defined E0 in a way that doesn't require rewriting code from the base pygemini or mat_gemini repo?

[jdiazpena]

Exactly, it is not hard, I did most of it, but we need a default function (or simply a function that gives a constant value) for when a function for E0 is not given. I did not do that. I attach here my versions of functions that define Q and E0

`import numpy as np
import xarray

def Qsaid(pg: xarray.Dataset, Qpeak: float, Qbackground: float) -> xarray.Dataset:
mlon_mean = pg.mlon.mean().item()
mlat_mean = pg.mlat.mean().item()

beta=0.15
T=1/42
f=pg.mlon.data[:, None] - mlon_mean
shapelon=0*f

for i in range(np.size(f,0)):
    if abs(f[i])<(1-beta)/(2*T):
        shapelon[i]=1
    elif (1-beta)/(2*T)<abs(f[i]) and abs(f[i])<(1+beta)/(2*T):
        shapelon[i]=0.5*(1+np.cos( (np.pi*T/beta)*(abs(f[i])-(1-beta)/(2*T)) ))
    else:
        shapelon[i]=0



if "mlat_sigma" in pg.attrs:
    shapelat = np.exp(
        -((pg.mlat.data[None, :] - mlat_mean - 8.65 * pg.mlat_sigma) ** 2) / (2 * pg.mlat_sigma**2)
    )
else:
    raise LookupError("precipation must be defined in latitude, longitude or both")

Q=Qpeak*shapelon*shapelat

Q[Q < Qbackground] = Qbackground

return Q

`

`import numpy as np
import xarray

def Esaid(pg: xarray.Dataset, Epeak: float, Ebackground: float) -> xarray.Dataset:
mlon_mean = pg.mlon.mean().item()
mlat_mean = pg.mlat.mean().item()

beta=0.15
T=1/42
f=pg.mlon.data[:, None] - mlon_mean
shapelon=0*f

for i in range(np.size(f,0)):
    if abs(f[i])<(1-beta)/(2*T):
        shapelon[i]=1
    elif (1-beta)/(2*T)<abs(f[i]) and abs(f[i])<(1+beta)/(2*T):
        shapelon[i]=0.5*(1+np.cos( (np.pi*T/beta)*(abs(f[i])-(1-beta)/(2*T)) ))
    else:
        shapelon[i]=0



if "mlat_sigma" in pg.attrs:
    shapelat = np.exp(
        -((pg.mlat.data[None, :] - mlat_mean - 8.65 * pg.mlat_sigma) ** 2) / (2 * pg.mlat_sigma**2)
    )
else:
    raise LookupError("precipation must be defined in latitude, longitude or both")

E=Epeak*shapelon*shapelat

E[E < Ebackground] = Ebackground

return E

`

[mattzett]

Awesome thanks; I'm going to leave this open until we get that extra code in place as a reminder that this functionality needs to be added.

[mattzett]

Just as a refresh; @hayleyclev and I have decided that there needs to be a means to have Qprecip_function to return both a total energy flux and a characteristic energy to accommodate a more general set of inputs.