% QDQ - Cumulative heat function and its derivative
%  [Q,DQ] = QDQ(THETA)
% THETA must be in radians. Function is correctly vectorized.
% Reads global variables THETA_S, THETA_B, QBAR

% Philip D Loewen, 2010-02-24

function [Q,dQ] = qdq(theta)

global global_THETA_S global_THETA_B global_QBAR
THETA_S = global_THETA_S;
THETA_B = global_THETA_B;
QBAR    = global_QBAR;

% Shift theta by 2n(pi) to put all values into [-pi,pi]
ThetaFixer = round(theta/2/pi/(1+eps));
Theta      = theta - 2*pi*ThetaFixer;

% Initialize return vector with all zeros
% (This is the correct value for small theta anyway.)
Q  = zeros(size(Theta));
dQ = zeros(size(Theta));

% Put global value QBAR into large-angle positions for Q
% (Note that Q'=0 also for these positions.)
Q(Theta>(THETA_S+THETA_B))  = QBAR;

% Construct a logical 0/1 vector that selects angles in the burn region
steep_bit = (Theta > THETA_S) & (Theta < (THETA_S+THETA_B));

% Apply fancy formulas in the burniregion
Q(steep_bit)  = 0.5 * QBAR * (1 - cos(pi*(Theta(steep_bit)-THETA_S)/THETA_B));
dQ(steep_bit) = 0.5 * QBAR * (pi/THETA_B) * sin(pi*(Theta(steep_bit)-THETA_S)/THETA_B);