s = poly(0,'s');

// Active feedback network
c15 = 3.3e-6;
c12 = 330e-9;
R7 = 180;

z1 = 1/s/c15 + R7;
z2 = 1/s/c12;

// Input RC (pole)
R4=1.8e3;
R5=1.5e3;
C13=47e-9;

// Output RC (pole)
R6=220;
C14=100e-9;

// Loop filter transfer function = f1*f2*f3
f1 = 1 / (R4+R5) / (1+s*C13/(1/R4+1/R5));
f2 = 1 / (1/z1+1/z2);
f3 = 1 / (1+s*R6*C14);

// Split 2nd-order term into two 1st-order terms
f2 = f2 * s;
f4 = 1/s;

// Convert to state space
f1ss = tf2ss(syslin('c', f1), 1e-40);
f2ss = tf2ss(syslin('c', f2), 1e-40);
f3ss = tf2ss(syslin('c', f3), 1e-40);
f4ss = tf2ss(syslin('c', f4), 1e-40);

format(20,'e');
lines(0);
clc;

Ts=1e-6;      // Timestep
PreWarp=500;  // 500 Hz pre-warping

// Continuous --> Discrete
dss1 = cls2dls(f1ss, Ts, PreWarp);
dss2 = cls2dls(f2ss, Ts, PreWarp);
dss3 = cls2dls(f3ss, Ts, PreWarp);
dss4 = cls2dls(f4ss, Ts, PreWarp);

// Back to transfer functions
dtf1 = ss2tf(dss1, 1e-40);
dtf2 = ss2tf(dss2, 1e-40);
dtf3 = ss2tf(dss3, 1e-40);
dtf4 = ss2tf(dss4, 1e-40);

// Zeros
n1 = coeff(numer(dtf1/dss1(5)));
n2 = coeff(numer(dtf2/dss2(5)));
n3 = coeff(numer(dtf3/dss3(5)));
n4 = coeff(numer(dtf4/dss4(5)));

// Poles
d1 = coeff(denom(dtf1/dss1(5)));
d2 = coeff(denom(dtf2/dss2(5)));
d3 = coeff(denom(dtf3/dss3(5)));
d4 = coeff(denom(dtf4/dss4(5)));

printf("POLEZERO pz1(%.30gL, %.30gL);\n", d1(1), n1(1));
printf("POLEZERO pz2(%.30gL, %.30gL);\n", d2(1), n2(1));
printf("POLEZERO pz3(%.30gL, %.30gL);\n", d3(1), n3(1));
printf("POLEZERO pz4(%.30gL, %.30gL);\n", d4(1), n4(1));
printf("%.30gL\n", dss1(5)*dss2(5)*dss3(5)*dss4(5));