You can solve it as follows:
Assuming a stable circuit in the Laplace domain we have Vn = Vp = GND.
$$
I_D - I_{CD} + I_{Rf} + I_{CF} + I_n = 0\\
I_{Rf}=\frac{V_{o2}-V_n}{R_f}\\
I_{Cf}=\frac{V_{o}-V_n}{\frac{1}{sC_f}}\\
V_{o2}=\frac{Z_L}{Z_L+R_{Fi}}Vo\\
Z_L=\frac{1}{sC_{Fi}}//R_f//R_L
$$
where Vo2 is the output after RFilter and ZL is the load impedance after RFilter. (I've denoted CFilter and RFilter as CFi and RFi)
Assuming In = 0 and neglecting CD, as it is not part of the filter circuit, and substituting it all into the node equation we can solve for Vo.
If we divide by ID (our input current), we get the solution for Vo:
$$
\frac{V_o}{I_D}=-\frac{\mathit{R_{Fi}}\cdot \mathit{R_f}+\left( \mathit{C_{Fi}}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot s+\mathit{R_f}+\mathit{R_{Fi}}\right) \cdot \mathit{R_L}}{\left( 1+\left( \mathit{C_f}\cdot \mathit{R_f}+\mathit{C_f}\cdot \mathit{R_{Fi}}\right) \cdot s+\mathit{C_{Fi}}\cdot \mathit{C_f}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot {{s}^{2}}\right) \cdot \mathit{R_L}+\mathit{C_f}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot s}
$$
With the known relation between Vo and Vo2 we can determine
$$
\frac{V_{o2}}{I_D}=-\frac{\mathit{R_f}\cdot \mathit{R_L}}{\left( 1+\left( \mathit{C_f}\cdot \mathit{R_f}+\mathit{C_f}\cdot \mathit{R_{Fi}}\right) \cdot s+\mathit{C_{Fi}}\cdot \mathit{C_f}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot {{s}^{2}}\right) \cdot \mathit{R_L}+\mathit{C_f}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot s}
$$
which is your second order filter equation.
To incorporate CD is a bit difficult because you get a current divider between your input impedance of the filter and CD.
But as you can see from the results above, it is not part of the second order filter.
Edit: I just realized that for an analytical calculation under the assumption of infinite open loop gain, CD is shorted by the input impedance, which approaches 0 (because Vn=Vp=0V=AC GND). Therefore it seems reasonable to neglect it. If you want to incorporate CD, you would have to use the finite (frequency dependent) open loop gain of the opamp to solve the circuit. Needless to say, that this makes the calculations even more lengthy.
Hope the answer helps you. You best check my result as the calculation is a bit tedious and prone to errors.
As for your first question: If you need a defined output impedance, you might want to add another unity gain amplifier as line driver.
Update (v2): (in the first version I still used the voltage divider for Vo2, that was wrong, since Vn is now unknown)
Your schematic as well as your calculations are (to my best knowledge) correct up to (and not including)
$$ V_{o2} = V_o\cdot \frac{Z_{fi}\|R_{f}}{R_{fi}+Z_{fi}\|R_{f}} $$
As detailed below, this formula is valid only for infinite open loop gain.
\$R_L\$ was indeed the coax load.
Then you use a voltage divider for \$V_{o2}\$, but don't take into account that Vo is a driven potential, directly related to the the input Vn (and Vp). In other words, your circuit is missing the (unknown) output impedance of the opamp. If you add it, the formula \$ V_{o2}= V_n\cdot \frac{Z_{fi}}{(Z_f+R_{fi})\|R_f+Z_{fi}}\$ doesn't match the circuit.
The way I would solve it is as follows:
Use the equations from my first solution (or yours), replace \$I_D\$ with an generic input current \$I_{in}\$, add the equation
$$V_o = a_f (V_p - V_n) = - a_f V_n$$
You also have to replace the equation for \$V_{o2}\$, since the current divider is now to Vn and not to Gnd. We solve that by adding another node equation:
$$I_{R_{Fi}}-I_{C_{Fi}}-I_{R_f}=0;$$
and branch equations (my directions):
$$I_{R_{Fi}}=\frac{(V_o-V_{o2})}{R_{Fi}}\\
I_{C_{Fi}}=\frac{V_{o2}}{\frac{1}{s CFi}}
$$
Substituting into the new node equation you can solve for Vo2:
$$
V_{o2}=\frac{\mathit{R_{Fi}}\cdot \mathit{V_n}+\mathit{R_f}\cdot \mathit{V_o}}{\mathit{C_{Fi}}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot s+\mathit{R_f}+\mathit{R_{Fi}}}
$$
With this you add one unknown (\$V_n\$) and one equation, therefore the next step is to solve for \$V_n\$, then use it to determine \$V_o\$, and then \$V_{o2}\$ (without any unknowns remaining).
You will retain \$a_f\$ as a variable throughout the calculation. You can check the result as \$
\mathop {\lim }\limits_{a_f \to \infty ; \; R_L \to \infty } V_{o2}\$ should give the previous result.
With this, you can calculate the input impedance by using $$Z_{in}=\frac{V_n}{I_{in}}$$
This will approach 0 for \$a_f \rightarrow \infty \$
From this you can calculate your input currents:
$$I_{in} = \frac{Z_{CD}}{Z_{CD}+Z_{in}} I_D\\
I_{CD} = \frac{Z_{in}}{Z_{CD}+Z_{in}} I_D
$$
If you now substitute the first equation for \$I_{in}\$ in your previous result, you get the exact analytical solution as you requested, incorporating all elements your circuit has.
Only downside is: now you need to model \$a_f\$. Maybe a first order low pass model will be sufficient. Substitute it into your solutions and you are done.
I would carefully check your results agains e.g. a spice simulation or real world results, as there are many non-idealities (input bias currents, offsets) that are still not part of this analysis.
You should get:
$$
\frac{V_{o2}}{I_{in}}=-\frac{\mathit{af}\cdot \mathit{R_f}- \mathit{R_{Fi}}}{\left( 1+\mathit{af}\right) \cdot \mathit{C_{Fi}}\cdot \mathit{C_f}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot {{s}^{2}}+\left( \left( \left( 1+\mathit{af}\right) \cdot \mathit{C_f}+\mathit{C_{Fi}}\right) \cdot \mathit{R_{Fi}}+\left( \mathit{af}+1\right) \cdot \mathit{C_f}\cdot \mathit{R_f}\right) \cdot s+\mathit{af}+1}
$$ and $$
Z_{in}=\frac{\mathit{R_{Fi}}+\mathit{R_f}+\mathit{C_{Fi}}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot s}{\left( 1+\mathit{af}\right) \cdot \mathit{C_{Fi}}\cdot \mathit{C_f}\cdot \mathit{R_{Fi}}\cdot \mathit{R_f}\cdot {{s}^{2}}+\left( \left( \left( 1+\mathit{af}\right) \cdot \mathit{C_f}+\mathit{C_{Fi}}\right) \cdot \mathit{R_{Fi}}+\left( \mathit{af}+1\right) \cdot \mathit{C_f}\cdot \mathit{R_f}\right) \cdot s+\mathit{af}+1}
$$
Good luck.
I believe the authors used fast analytical circuits techniques or FACTs but I am not sure about their results. The principle is based on the generalized Extra-Element Theorem or EET from Dr. Middlebrook (https://en.wikipedia.org/wiki/Extra_element_theorem). If I try to determine the transfer function of the below passive circuit in which I put arbitrary components values:
I can write the transfer function linking B to A under the form \$H(s)=H_0\frac{N(s)}{D(s)}\$ in which \$H_0\$ is the gain determined for \$s=0\$ when all caps are open. The first thing is to open the caps and determine the transfer function in this mode. Then, reduce the excitation to 0 V (replace it by a short circuit) and "look" at the resistance from the capacitor terminals in this mode. This resistance multiplied by the capacitor forms the time constant, \$\tau=RC\$. Summing these time constants gives \$b_1\$ in \$D(s)\$. \$b_2\$ is obtained by summing a product of time constants re-using one of the time constants in \$b_1\$. If I chose to re-use \$\tau_L\$ meaning the time constant associated with \$C_L\$, then I will put this capacitor in its hi-frequency state (a short circuit) and "look" at the resistance offered by \$C_f\$ terminals in this configuration. All these operations are shown in the below sketch:
From there, you can assemble a second-order polynomial form:
\$D(s)=1+s(\tau_1+\tau_2)+s^2\tau_1\tau_{12}=1+\frac{s}{\omega_0Q}+(\frac{s}{\omega_0})^2\$ calculate the quality factor \$Q\$ and factor \$D\$ as two cascaded poles if \$Q<<1\$: \$D(s)=(1+\frac{s}{\omega_{p1}})(1+\frac{s}{\omega_{p2}})\$ in which \$\omega_{p1}=Q\omega_0\$ and \$\omega_{p2}=\frac{\omega_0}{Q}\$.
For the zeros, you can either use a null double injection (NDI) or a generalized 2nd-order form as described in https://www.amazon.com/Linear-Circuit-Transfer-Functions-Introduction/dp/1119236371/ref=asap_bc?ie=UTF8. It is a bit longer than the NDI but sometimes simpler to implement. The corresponding schematic for this exercise is here:
If you develop the numerator, you obtain a double zero:
\$N(s)=1+\frac{s}{\omega_{0N}Q_N}+(\frac{s}{\omega_{0N}})^2\$ and cannot really split it as two cascaded zeros as the two are coincident (\$Q\$ is almost 1)
The final transfer function is given in the below Mathcad shots with the frequency response. You can see that this circuit builds phase boost between the poles and zeros and will certainly improve phase margin at crossover when used in a compensator.
The FACTs are truly valuable in this case because you can determine the transfer function by inspection, without writing a line of algebra. And if you make a typo, you can correct the individual sketch that is causing problems. You have a tutorial taught at APEC in 2016 that you can use for a smooth introduction to the technique: http://cbasso.pagesperso-orange.fr/Downloads/PPTs/Chris%20Basso%20APEC%20seminar%202016.pdf.
Best Answer
Imagine a series resistance Rs in S1 normally closed contact. The current If flowing through the feedback resistor Rf1 will result in the voltage at the op-amp output being higher than the ideal voltage by Rs*If.
But the voltage at the right hand side of Rf1 will be unaffected (because it's inside the feedback loop).
If you add the second switch you can pick off that voltage and provided there is negligible current through the resistance of S2 (to the load) you will have eliminated all the error caused by the switch resistance.
For illustration, consider the below example where you have two switches with different resistances in each position.
simulate this circuit – Schematic created using CircuitLab
This dual-switch configuration is frequently useful when you have high resistance switches such as CMOS analog switches.