% 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);