Error propagation when dividing by fitted model

#1
Hoping someone can help me / point me in the right direction with this one. Seems like it should be straight forward, but just can't quite get my head around it or find the info I need.

BACKGROUND:

I am looking at an astronomical radio data set. For a bunch of radio sources, I have measurements of the total radio intensity in ~100 different frequency channels, and the intensities of the orthogonal linear polarization in the same frequency channels. I also have associated measurements of the noise level in each channel, which is Gaussian.

Since I am interested in the frequency-dependent fractional polarization, I need to divide the linearly polarized intensities by the total intensity. Rather than do this on a channel by channel basis, I divide the linearly polarized intensities in each frequency channel by the value of a model fitted to the total intensity data across the full frequency range, evaluated at that frequency. The fitted model is a power law - i.e. a straight line fitted to the log(total intensity) against log(frequency), and represents a good fit to the data.

QUESTION:

My question is this - I have a bunch of data points with associated errors (the linearly polarized intensities), and I have divided these by a model (a fit to the total intensities) which has errors associated with its slope and y-intercept. How do I do the error propagation in this case? In other words, given that I'm dividing a set of discrete data points (with their associated errors) by a model fitted to another set of discrete data points (with errors associated with its slope and y-intercept) to get a final set of data points, how do I propagate the errors through to find the resultant error on these final data points?