Yes, you can get rid of both loops (but you can't use powerlaw_fit). You can exploit the fact that the log of your power law is a linear function of log(wavelength). In other words:
y = (wavelength/a)^b
becomes
log(y) = b*log(wavelength) + constant
Then, you could solve for your parameters using matrix operations.
To simplify the problem a bit, let's pretend that your nl by ns detector array is instead a 1D detector with N = nl * ns pixels. Each pixel, i, contains your 9 measured data values as a function of wavelength for that pixel, y_i = [ y_i_1 y_i_2 ... y_i_9 ].
Now, we take the log of both sides of the power law function. Let y'_i = log(y_i), wl' = log(wavelength), and a'_i = -b_i * log(a_i), where a_i and b_i are the power law parameters we are trying to find. We can then encapsulate everything into a matrix equation:
[ 1 wl'_1 ] [ a'_1 a'_2 ..... a'_N ] [ y'_1_1 y'_2_1 ... y'_N_1 ]
[ 1 wl'_2 ] [ b_1 b_2 ..... b_N ] = [ y'_1_2 y'_2_2 ... y'_N_2 ]
[ 1 wl'_3 ] [ y'_1_3 y'_2_3 ... y'_N_3 ]
... ...
[ 1 wl'_9 ] [ y'_1_9 y'_2_9 ... y'_N_9 ]
That first matrix with the wavelengths is the design matrix, M, and let's call the parameter matrix P, and the data matrix, Y. To restate the above more compactly,
MP = Y
To solve for P, we multiply both sides by the transpose of M (denoted by M_T), and then multiply both sides of the resulting equation by the inverse of M_T M (denoted by (M_T M)^-1). That gives:
P = (M_T M)^-1 M_T Y
This doesn't incorporate errors, but I'm pretty sure there's a way to do that. One note about errors is that doing this log-trick does mess with the distribution of the errors, so either take care to do that part properly, or pretend I didn't say anything!
Anyway, I wrote this code to generate fake data and then recover the parameters using the linear algebra described above. In the case of this noiseless data, at least, the method works and it's very fast. I tried to comment the code sufficiently, but let me know if you need any clarification. Here is it:
function simple_powerlaw, wl, a, b
; Returns y values for power law
; WL - your wavelength vector
; A - amplitude coefficient
; B - exponent coefficient
return, (wl/a)^b
end
------------------------------------------------------------
function make_fake_data, real_a, real_b
; Returns an Npixel x Nwl array, where each column is the
; measured data as a function of wavelength for a given pixel
; define dimensions of detector
ncol = 100L
nrow = 150L
npix_total = ncol * nrow
; your wavelength data
wavelength = [412,442,488,530,554,645,747,857,1242]
nwl = n_elements(wavelength)
; Initialize array of data (Npixels x Nwl)
fake_data = dblarr(npix_total, nwl)
; Generate the power-law data at each pixel
for pix = 0, npix_total-1 do begin
; get the power law coefficients for the given pixel
a_i = real_a[pix]
b_i = real_b[pix]
; generate the power-law y-values as a function of wavelength
; for the given pixel
fake_data[pix, *] = simple_powerlaw(wavelength, a_i, b_i)
endfor
return, fake_data
end
-----------------------------------------------------------
pro plaw_to_mlr
;;; MAKE FAKE DATA
; define the numbers of rows and columns
ncol = 100
nrow = 150
npix_total = ncol * nrow
; wavelength from your post
wavelength = [412,442,488,530,554,645,747,857,1242]
nwl = n_elements(wavelength)
; Randomly generate power-law parameters
; Each pixel will have two parameters, defined by the following vectors
; I arbitrarily define the amplitudes to range from 0-5, and the
; exponents to range from 1-2, all randomly assigned
; Also, since we are making a fit at each pixel, I'm treating the 2D
; detector as a 1D array of pixels, for convenience
real_a = randomu(seed, npix_total) * 5d0
real_b = randomu(seed, npix_total) + 1d0
; make the fake data
data = make_fake_data(real_a, real_b)
;;; DO LINEAR ALGEBRA
; [design_matrix] x [parameter_matrix] = [data_matrix]
; Solve to get => [parameter_matrix] = ([design_matrix]^T x [design_matrix])^-1 x [data_matrix]
; Design matrix (a column of 1's and a column of the (log) wavelength vector
logwavelength = alog10(wavelength)
design_matrix = [transpose(replicate(1d0, nwl)), transpose(logwavelength)]
; y-value/data matrix
logdata_arr = alog10(data)
; get transpose of design matrix (dm), and dmT x dm
dmT = transpose(design_matrix)
dmTdm = dmT ## design_matrix
; get inverse of dmTdm
dmTdm_inv = invert(dmTdm)
; Solve for the parameters
; This array is a Npix x 2 array
; The second row is an Npix vector corresponding to the value of b at
; each pixel
; The first row contains a constant term, which is related to a: const = -b * log a
fit_parameters = (dmTdm_inv ## dmT) ## logdata_arr
fit_b = fit_parameters[*,1]
fit_a = 10d0^(-1d0*fit_parameters[*,0]/fit_b)
; Check against 'true' parameter values
print, total(real_a - fit_a)
print, total(real_b - fit_b)
end
As it is written, the code just prints sum total of the differences between the 15,000 true parameter values and the 15,000 derived parameter values, for a and b, respectively. Here's the output of the call:
IDL> plaw_to_mlr
% Compiled module: PLAW_TO_MLR
0.028222162
0.0032718051
