However, I need p[4] to be strictly between 0 and pi. Is it possible to do it?
In addition, is there a wat to calculate or plot the 95% confidence intervals of the fit?

If not, is there another recommended package that can deal with this effectively?

I don’t think LsqFit.jl supports bounded optimization. However in this case it would not be too hard to transform the parameters found by the optimization to the desired interval.