function m=zhadmtrx(tempk, tempq, tempv, theta) % Transition matrix for BPS shaker model. Model parameters are: rate % parameters "tempk", charge parameters "tempq". "tempv" is the % applied voltage temperature=12; % Temperature in Celsius degrees kT=(273.15+temperature)*0.08625; A = tempk(1)*exp(tempq(1)*tempv/kT); B = tempk(2)*exp(tempq(2)*tempv/kT); g = tempk(3)*exp(tempq(3)*tempv/kT); d = tempk(4)*exp(tempq(4)*tempv/kT); OC = tempk(5)*exp(tempq(5)*tempv/kT); CO = tempk(6)*exp(tempq(6)*tempv/kT); thet=theta; m =[-4*A, B, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 4*A, -B-3*A-g, 2*B, 0, 0, d, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 3*A, -2*B-2*A-2*g, 3*B, 0, 0, d, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 2*A, -3*B-A-3*g, 4*B, 0, 0, d, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, A, -4*B-4*g, 0, 0, 0, d, 0, 0, 0, 0, 0, 0, 0; 0, g, 0, 0, 0, -d-3*A, B, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 2*g, 0, 0, 3*A, -B-d-2*A-g, 2*B, 0, 2*d, 0, 0, 0, 0, 0, 0; 0, 0, 0, 3*g, 0, 0, 2*A, -2*B-d-A-2*g, 3*B, 0, 2*d, 0, 0, 0, 0, 0; 0, 0, 0, 0, 4*g, 0, 0, A, -3*B-d-3*g, 0, 0, 2*d, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, g, 0, 0, -2*d-2*A, B, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 2*g, 0, 2*A, -B-2*d-A-g, 2*B, 3*d, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 3*g, 0, A, -2*B-2*d-2*g, 0, 3*d, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, g, 0, -3*d-A, B, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2*g, A, -B-3*d-g, (4*d)/thet, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, g, -((4*d)/thet)-OC, CO; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, OC, -CO];