I am trying to fit some data points $d_i$ to a non-linear model function $m_i$, which depends on a number of fit parameters $f_k$ (I want to determine these) and also on some known, constant values $c_j$. For this I am using the Levenberg–Marquardt (LM) algorithm, which iteratively minimizes
\begin{equation} \chi^2=\sum_i\left[\frac{d_i-m_i\left(f_1,f_2,\ldots,c_1,c_2,\ldots\right)}{\sigma_i\left(d_i\right)}\right]^2\,,\tag{1} \end{equation}
where $\sigma_i^2\left(d_i\right)$ is the variance of $d_i$ (I'll call $\sigma_i$ the error of $d_i$). I have an estimate for the errors, so let's say they are known.
Now, I want to deal with the "constant" values $c_j$, which actually come from different measurements and thus have some (known) error themselves. My first idea is to adjust the term in the denominator, with the simplest type of error propagation:
\begin{equation} \chi^2=\sum_i\left[\frac{d_i-m_i\left(f_1,f_2,\ldots,c_1,c_2,\ldots\right)}{\sqrt{\sigma_i^2\left(d_i\right)+\sum_j\left(\frac{\partial{}m_i}{\partial{}c_j}\sigma\left(c_j\right)\right)^2}}\right]^2\,.\tag{2} \end{equation}
My question is how to calculate the partial derivatives $\frac{\partial{}m_i}{\partial{}c_j}$ correctly, since they depend on the current value of the fit parameters $f_k$. I see two solutions, but I don't find either very satisfying:
1) Apply the LM algorithm multiples times. For example: Calculate the derivatives with a starting guess $f_k^0$ and execute the LM fit, giving the solution $f_k^1$. Then, use $f_k^1$ to calculate the derivatives and execute a second LM fit.
2) Adjust the LM algorithm code to actually recalculate the denominator after each LM-iteration with the current best estimates of $f_k$.
How would you recommend I approach this problem? Any hints in general or reading suggestions are also appreciated.