analc_endsversion/AnalyticLC/MALookupTable.m
2022-04-11 15:28:22 +02:00

105 lines
2.9 KiB
Matlab

function Output = MALookupTable(varargin)
% MALookupTable: generate a look-up table for Mandel-Agol model values, or
% use the table and read from it
%Inputs: 'init' - calculate a table of the MA values for the given rVec,
%zVec, u1, u2
% 'load' = return all the table and the row/col values in a
% structure
% r,z = return the Mandel-Agol model values for the given r,z.
% u1 and u2 which will be used are the ones that were used for
% generating the table. z can be a vector, r must be a scalar.
% Yair Judkovsky, 25.10.2020
persistent rVec zVec u1 u2 R Z MATable
%return a specific Mandel-Agol model value
if isnumeric(varargin{1})
%get r, z values from the input
r = varargin{1};
z = varargin{2};
Requested_u1 = varargin{3};
Requested_u2 = varargin{4};
%if no table exists, or if the requested u1,u2 do not match the table, just calculate the values
if isempty(MATable) || Requested_u1~=u1 || Requested_u2~=u2
Output = occultquadVec(r,z,Requested_u1,Requested_u2);
return
end
Output = ones(size(z));
%old method - plug in "1" wherever outside the planet radius range
% Output(z<1+r) = interp2(R,Z,MATable,z(z<1+r),r,'linear',1);
%new method - outside the radius range calculate using occultquadVec. This
%is useful because now the radius range can be reduce, yielding a fast
%calculation for most of the points
ind = r>min(rVec) & r<max(rVec);
Output(z<1+r & ind) = interp2(R,Z,MATable,z(z<1+r & ind),r,'linear');
Output(z<1+r & ~ind) = occultquadVec(r,z(z<1+r & ~ind),u1,u2);
return
end
%perform the calculation. Table will be made such that each row represents
%a value of r, and each column represents a value of z
if isequal(varargin{1},'init')
%set range of r, z, u1, u2. Can also be read from varargin if
%exists, where the 2,3,4,5 arguments will be rVec, zVec, u1, u2.
rVec = 0.5e-2:2e-3:0.2;
zVec = 0:2e-3:1.25;
u1 = 0.36;
u2 = 0.28;
%read inputs from varargin, if they exist
if length(varargin)==5
rVec = varargin{2};
zVec = varargin{3};
u1 = varargin{4};
u2 = varargin{5};
end
%generate a table of r and z values
[R,Z] = meshgrid(zVec,rVec);
%set a table of model values for r and z
MATable = zeros(length(rVec),length(zVec));
for j = 1:length(rVec)
MATable(j,:) = occultquadVec(rVec(j),zVec,u1,u2);
end
fprintf('calculated Mandel-Agol values \n');
Output = [];
return
end
%load the table
if isequal(varargin{1},'load')
Output.rVec = rVec;
Output.zVec = zVec;
Output.u1 = u1;
Output.u2 = u2;
Output.R = R;
Output.Z = Z;
Output.MATable = MATable;
return
end
%clear the values
if isequal(varargin{1},'clear')
clear rVec zVec u1 u2 R Z MATable
Output = [];
return
end