I have an equation for the reflection coefficient of a certain circuit (I can give more context if you are interested), and it is given by
$$\Gamma(\omega) = \frac{8 J_0^2 (-2 i \gamma +i \kappa +4 J_1+4 \omega -2 (\text{$\omega_0$}+\text{$\omega_1$}))+(i \gamma +2 J_1-2 \omega +2 \text{$\omega_1$}) \left(-4 J_1^2-(\gamma +2 i (\omega -\text{$\omega_1$})) (\gamma -\kappa +2 i (\omega -\text{$\omega_0$}))\right)}{8 J_0^2 (-2 i \gamma -i \kappa +4 J_1+4 \omega -2 (\text{$\omega_0$}+\text{$\omega_1$}))+(i \gamma +2 J_1-2 \omega +2 \text{$\omega_1$}) \left(-4 J_1^2-(\gamma +2 i (\omega -\text{$\omega_1$})) (\gamma +\kappa +2 i (\omega -\text{$\omega_0$}))\right)}$$
Definitely not a very pretty expression, but on the other hand not too complicated; the numerator and denominator only differ by a few signs. I should say that $J_0,J_1,\omega_0,\omega_1,\gamma,\kappa$ are all parameters, and that the variable is $\omega$. They are always real and positive, $J_0 > J_1$ and $\omega_0 \leq \omega_1$. This expression is complex, but when looking at either the absolute value or the argument, you can see three resonances, as illustrated in this figure

These figures were generated in Mathematica, using $\omega_0 = 47.0, \omega_1 = 47.3, J_0 = 0.56, J_1 = 0.026, \kappa = 0.00011, \gamma = 0.19$. This is not too relevant to my question, but in case you were interested.
Now, I have done an experiment on a system that exhibits dynamics consistent with this model, and I'm trying to fit it. I have the full complex dataset, so I can look at the absolute value, the argument, the polar part, whatever I think is easiest. My problem comes from the fact that there are a lot of parameters to be fit, and that this requires very, very good initial parameters to have any chance. So my question concerns this.
What I want to find out is if, given the equation above, I can think of an analytical way of finding starting values for the fitting routine, for the parameters. My first approach has been quite simple, as I have a lot of trouble actually extracting information from the equation I gave. I look at the pictures, and I see that the middle resonance is pretty much at $\omega_{0,1}$. So this is my starting value for these two. I then looked at the spacing between the first and second resonance, and see that this appears to be approximately $2J_0$. So that was my initial guess for that one. Then, I (again from how the experiment is set up) know that $J_1$ will usually be about a 10th of $J_0$. You can already see from the data I put in that this is not correct at all, but simply an order of magnitude.
But then the problem of $\kappa$ and $\gamma$ arises. To be honest, I am not sure what to do with these. I feel like they are related to the width and height of the peaks, but I can't see how. Perhaps this is something you can help me with. This is really the main part of the question; how can I interpret the role of $kappa$ and $gamma$, so that I can extract some basic guess from, say, the resonances in the data?
Moreover, my initial guess for the others can probably also be improved upon by actually understanding the equation. So yeah, my question boils down to if someone can help me understand the contribution of the terms in my equation, so that I can get better estimates of their values from the data.
Another thought I had: perhaps you have an opinion on what would be better to use for a fit. Absolute value, argument, polar plot?
Oh, I'd like to add, I'm not sure what kind of tag to put this under, so suggestions are welcome.
Additional information
I cannot judge how important it is, but I can talk about where the equation comes from. Feel free to disregard most of this if you think it does not help make the question easier. My question is mainly about understanding the equation given above, but some context might help.
We have a system of four LC resonators that are coupled in a periodic fashion. They each have a resonant frequency of their own (if you don't look at the system as a whole), which is $\omega_i$. One of them is connected to the outside world, hence its resonance frequency $\omega_0$ is different from the other three, which get $\omega_1$. Now, the coupling strength between nearest neighbours is $J_0$, and between next-nearest neighbours is $J_1$.
This system has a Hamiltonian
$$\left( \begin{array}{cccc} \text{$\omega_0$} & J_0 & J_1 & J_0 \\ J_0 & \text{$\omega_1$} & J_0 & J_1\\ J_1 & J_0 & \text{$\omega_1$} & J_0 \\ J_0 & J_1 & J_0& \text{$\omega_1$} \\ \end{array} \right)$$
The eigenvalues and eigenvectors are rather complicated. If we instead set $\omega_0=\omega_1$ (they will often be rather close as the asymmetry comes from coupling to the outside world, which is taken to be rather weak) we get simple eigenvalues: $\omega_1-J_1,\omega_1-J_1,\omega_1-2J_0+J_1,\omega_1+2J_0+J_1$. Now that I write this out, these values do seem to correspond to the locations of the peaks, so I can probably do some elimination here and get better starting values for these parameters.
It becomes more complicated when one considers where the reflection coefficient comes from. Without going into the physics, it comes from solving a system $Ma = b$, where $M$ is the matrix
$$M = \left( \begin{array}{cccc} \frac{i \gamma }{2}+\frac{i \kappa }{2}-\omega +\text{$\omega_0$} & J_0 & J_1 & J_0 \\ J_0 & \frac{i \gamma }{2}-\omega +\text{$\omega_1$} & J_0 & J_1 \\ J_1 & J_0 & \frac{i \gamma }{2}-\omega +\text{$\omega_1$} & J_0\\ J_0 & J_1 & J_0 & \frac{i \gamma }{2}-\omega +\text{$\omega_1$} \\ \end{array} \right)$$
and $b$ the vector $$b = \left\{i a_{in} \sqrt{\kappa },0,0,0\right\}$$
One then computes the equation I gave through $$\Gamma(\omega) = 1-\left\{\sqrt{\kappa },0,0,0\right\}a/a_{in}$$
What exactly $a_{in}$ is is irrelevant, it always cancels out in the end. So this gives you the final expression, and perhaps somewhere in here there is information on how to determine initial values for $\gamma,\kappa$, but I personally think this is a lot less probable than with the Hamiltonian, and their initial values will probably have to be estimated from the equation I gave at the beginning.
Oh, and as requested, the code used for the above figures is
GammaFunClosedKappa[\[Omega]_?NumericQ, \[Gamma]_, \[Kappa]_, J0_,
J1_, \[Omega]0_, \[Omega]1_] := (-(4 J1^2 + (\[Gamma] - \[Kappa] +
2 I (\[Omega] - \[Omega]0)) (\[Gamma] +
2 I (\[Omega] - \[Omega]1))) (2 J1 + I \[Gamma] -
2 \[Omega] + 2 \[Omega]1) +
8 J0^2 (4 J1 - 2 I \[Gamma] + I \[Kappa] + 4 \[Omega] -
2 (\[Omega]0 + \[Omega]1)))/(-(4 J1^2 + (\[Gamma] + \[Kappa] +
2 I (\[Omega] - \[Omega]0)) (\[Gamma] +
2 I (\[Omega] - \[Omega]1))) (2 J1 + I \[Gamma] -
2 \[Omega] + 2 \[Omega]1) +
8 J0^2 (4 J1 - 2 I \[Gamma] - I \[Kappa] + 4 \[Omega] -
2 (\[Omega]0 + \[Omega]1)))