function [P, Ec_int, Ei_int, Ec_gpts, Ei_gpts, nu_one, nu_two] = ... polar_func(params, E, Delta_t, freq_dep_yes, density_type) %%% Set the number of field points and the number of quadrature intervals. Nk = length(E); Np = 20; Nq = 20; Ni = 4*Np; Nj = 4*Nq; %%% Extract the parameters. %%% In the case of general densities, the integration limits are fixed %%% parameters. Otherwise they are determined dynamically from the %%% lognormal/normal or normal/normal parameters. if density_type == 3 Ec_int = params(1); Ei_int = params(2); params = params(3:end); end; if freq_dep_yes == 1 jump_freq = params(1); % jump_freq = 1/Tau beta_param = params(2); % beta_param = sqrt(eta*V/k*T) params = params(3:end); end; eta = params(1); P_R = 1; if density_type == 1 C = params(2); b_one = params(3); c_one = sqrt(C); b_two = params(4); c_two = sqrt(C); Ec_int = sqrt(b_one*7); Ei_int = sqrt(b_two*7); elseif density_type == 2 C = params(2); Ec_bar = params(3); c = params(4); c_one = sqrt(C); b_two = params(5); c_two = sqrt(C); Ec_int = Ec_bar*200^c; Ei_int = sqrt(b_two*7); end %%% Setup the quadrature points Ec_gpts and Ee_gpts along %%% with the densities nu_one and nu_two. hp = Ec_int/Np; v = gauss_weights(Np,hp); Ec_gpts = gauss_points(Np,hp,0); if density_type == 1 nu_one = c_one*exp(-((Ec_gpts).^2)/b_one); elseif density_type == 2 nu_one = c_one*exp(-(log(Ec_gpts/Ec_bar)/(2*c)).^2); elseif density_type == 3 nu_one = params(2:2+Ni-1); end hq = 2*Ei_int/Nq; % factor of two b/c this is neg-inf to pos-inf integrand w = gauss_weights(Nq,hq); Ei_gpts = gauss_points(Nq,hq,0-Ei_int); if density_type == 1 | density_type == 2 nu_two_pos = c_two*exp(-((Ei_gpts(Nj/2+1:Nj)).^2)/b_two); elseif density_type == 3 nu_two_pos = params(2+Ni:2+Ni+Nj/2-1); end nu_two_neg = rot90(rot90(nu_two_pos)); nu_two = [nu_two_neg; nu_two_pos]; %%% Setup the initial vectors V and W and the matrices epsilon_c and epsilon_k V = v.*nu_one; W = w.*nu_two; eps_c = Ec_gpts*ones(1,Nj); eps_k = ones(Ni,1)*Ei_gpts'; %%% Determine the polarization. P = zeros(Nk,1); if freq_dep_yes == 0 Delta = ones(Ni,Nj); Delta(:,1:Nj/2) = -1; P_bar = (1/eta)*eps_k + P_R*Delta; P(1) = V'*P_bar*W; %%% Start iteration for k=2:Nk dE = E(k) - E(k-1); eps_k = eps_k + dE; Delta = sign(eps_k + eps_c.*Delta); P_bar = (1/eta)*eps_k + P_R*Delta; P(k) = V'*P_bar*W; end; elseif freq_dep_yes == 1 P_I = P_R - eps_c/eta; x_pos_prev = ones(Ni,Nj); x_pos_prev(:,1:Nj/2) = 0; epsilon = 1e-3; fudge = 1e-50; % Start iteration for k = 2:Nk E_cur = E(k) + eps_k; P_cur_min_L = E_cur/eta - P_R; P_cur_min_R = E_cur/eta + P_R; r_neg_pos = (sign(eps_c - E_cur) + 1)/2.*... (erfc((P_cur_min_L - (-P_I + epsilon))/(beta_param*sqrt(2))) -... erfc((P_cur_min_L - (-P_I - epsilon))/(beta_param*sqrt(2))))./... (erfc((P_cur_min_L - (-P_I + epsilon))/(beta_param*sqrt(2))) + fudge) + (1 - sign(eps_c - E_cur))/2; r_pos_neg = (sign(eps_c + E_cur) + 1)/2.*... (erfc((P_I - epsilon - P_cur_min_R)/(beta_param*sqrt(2))) -... erfc((P_I + epsilon - P_cur_min_R)/(beta_param*sqrt(2))))./... (erfc((P_I - epsilon - P_cur_min_R)/(beta_param*sqrt(2))) + fudge) + (1 - sign(eps_c + E_cur))/2; p_neg_pos = jump_freq*r_neg_pos; p_pos_neg = jump_freq*r_pos_neg; x_pos_cur = (x_pos_prev + Delta_t*p_neg_pos)./(1 + Delta_t*(p_pos_neg + p_neg_pos)); x_pos_prev = x_pos_cur; P_neg_avg = -beta_param/sqrt(2*pi)*exp(-(-P_I - epsilon - P_cur_min_L).^2/(2*beta_param^2))./... (0.5*erfc(-(-P_I - epsilon - P_cur_min_L)/(beta_param*sqrt(2))) + fudge) + P_cur_min_L; P_pos_avg = beta_param/sqrt(2*pi)*exp(-(P_I + epsilon -P_cur_min_R).^2/(2*beta_param^2))./... (0.5*erfc((P_I + epsilon - P_cur_min_R)/(beta_param*sqrt(2))) + fudge) + P_cur_min_R; P_bar = x_pos_cur.*P_pos_avg + (1 - x_pos_cur).*P_neg_avg; P(k) = V'*P_bar*W; end; end;