function to generate calibration matrix for IQU data
0001 % function to generate calibration matrix for IQU data 0002 function calmat = polcal(data, model) 0003 0004 % data are corresponding values of I, Q, U, phi (parallactic angle) 0005 % as an n-by-4 matrix 0006 % model is a vector, [I (intensity) p (polarization fraction) theta (position angle)] 0007 % such that the model Stokes parameters are given by 0008 % I = I 0009 % Q = I * p * cos(2(theta + phi)) 0010 % U = I * p * sin(2(theta + phi)) 0011 % calmat is the 3-by-3 Mueller matrix that is the *inverse* of the 0012 % telescope response, ie if the true sky Stokes vector is S and the 0013 % observed vector is S', S = CS' 0014 % so what we want to minimise is 0015 % chisq = sum(mod-squared(CS' - S)) 0016 % as a function of C 0017 0018 calmat0 = eye(3); %starting value 0019 0020 % call fminunc using the nested function method 0021 [calmat, fval] = fminunc(@chisq, calmat0); 0022 0023 function f = chisq(calmat) 0024 f = 0; 0025 npoints = length(data(:,1)); 0026 for i = 1:npoints; 0027 phi = data(i,4); 0028 modelvec = [model(1), model(1)*model(2)*cosd(2*(model(3)+phi)), model(1)*model(2)*sind(2*(model(3)+phi))]; 0029 datavec = [data(i,1), data(i,2), data(i,3)]; 0030 chi = datavec*calmat - modelvec; 0031 f = f + norm(chi)^2; 0032 end 0033 end 0034 end 0035 %return calmat