I am currently working on a topic involving nuclear theory (the study of deformed nuclei). So for the past two days, I've been trying to fit some data; I developed a theoretical model that aims at describing the energy spectrum of a nucleus.

The analytical expressions for the energies have two free parameters, so I need to find those parameters such that the theoretical energies fit the experiment as much as possible.

I am used to Wolfram Mathematica, so I decided to solve this problem with the help of it.

The structure of my problem is like this: I have four energy levels (pairs of x,y values, where x is the spin and y is the energy). The model has an analytical expression for the excitation energy, and it depends on two free parameters and some other values which I know beforehand.

Now, all I need is to find the two free parameters and then represent the results with a series of plots.

I am using the `NonLinearModelFit`

function since the formulas are quite complex (and non-linear). With the obtained parameters, I go and give numerical values for the analytical expressions for the energy, the calculate the root mean square deviation (RMS) by using the well-known formula
$E_\text{RMS}=\sqrt{\frac{\sum_i(E_{exp}^i-E_{th}^i)^2}{N}}$.

After some effort, I've managed to get around 0.5 MeV. However, in our problem, such a deviation is quite poor. We aim at obtaining an RMS between 0.2-0.3 (0.35 at worst).
From what I can see, the issue is band-4, where half of the theoretical values are much smaller than the experimental ones.
I was also thinking about introducing some conditions in the fitting function (e.g. I already know that the second parameter `s`

must be positive) but I don't know which kind of conditions would help me here. You can also see that the analytical expressions for energies have an `If`

statement to them:

that is because when the `NonLinearModelFit`

starts to compute the best-fit parameters, it goes to some values for which the energies are complex. So to avoid that, I just make the energy a big number whenever its Imaginary part is different from zero. In this way, I can avoid the situation where the program runs for "useless" regions of the free parameters.

**Q**: Is there any way of manipulating that Fit function to decrease that RMS? (no more than 0.3) Or maybe just consider another approach rather than the current one?

Sorry for the long text, but I wanted to give you some context. The attached document contains all my code, with some explanations. The analytic expressions for the energies and everything else are irrelevant. The only essential part is the Fit function, which is marked by the green background.

Thank you in advance!

**Attachments:**