Re-express as a linear system Ac = b where c=[c_1, c_2] is a vector of the unknown coefficients, A is a 1000 \times 2 matrix whose columns are v_{app} u_s and \sin(\psi) \cos(\beta) /v_{app}, and b is the vector psi_dot.

Then solve with c = A\b. That’ll get you a least-squares solution to the problem, the one that minimizes sum res^2, using QR decomposition.

Absolutely not, performing linear regressions with LsqFit.jl is perfectly fine, and it is probably when the estimated standard deviations are the most accurate.

I think we are actually talking about a standard error of a parameter estimate here, not a standard deviation. You first should think about conceptually what you are actually trying to estimate here, as the docs of LsqFit say there are some assumptions being made

Assume the errors in each sample are independent, normal distributed with zero mean and same variance, i.e. ϵ∼N(0,σ2I)

which might or might not make sense (especially homoskedasticity is often a problematic assumption on real world data at least in the social sciences, I’m not a kite power plant modelling expert so YMMV).

But indeed your model is linear so you can just estimate it with GLM or FixedEffectModels and get heteroskedasiticty-robust standard errors as well if that matters to you.