The 2N3904 should work fine.
Assuming your 2 2N3904's match well (you can measure a number and select the 2 that match the best...):
At the output current peak, the output current should be related to the input current by: Iout = Iin/e
So, for Iout = 1mA, Iin = 1mA*e ~ 2.72mA
For Vin=12V, T=300K, estimate Vbe ~ .65V, R1 = (12V-0.65V)/2.72mA = 4.17k, so try 3.9k
At 300K, the thermal voltage is ~25.85mV, so RL = 25.85mV/2.72mA = 9.5, so try 10
I tried this for kicks with LTspice and got .96mA for Iout.
If the explanation in the link didn't work for you, I can try a slightly different approach.
For this approach, you need to feel comfortable with the Ebers-Moll model relating emitter current to Vbe and the differentiated version of that to create a relationship for gm. See Wikipedia (
https://en.wikipedia.org/wiki/Bipolar_junction_transistor#Ebers%E2%80%93Moll_model), if needed.
What it boils down to is that you don't want the output current to change as the input current changes. Ideally, Vbe2 doesn't change at all, but that's not possible with such a simple circuit. But what can we do?
Look at the equation for Vbe2:
Vbe2 = Vbe1 - Ic1*RL
How does Vbe2 change as Ic1 changes from a bias point, say due to a change in Vin?
For a bipolar transistor, there is a relationship between a change in Vbe and a change in Ic. This comes from the differentiation of the Ebers-Moll model:
gm = (change in Ic)/(change in Vbe) = Ic/VT (Ic == collector current; VT == thermal voltage)
gm, the transconductance, is a small signal parameter, so as the instantaneous Ic (total current) changes, gm will change, too; ignore that for now...
We'll use this gm relationship to help us fing RL:
What we want is to achieve a change in Vbe2, delta_Vbe2 = 0. Let's look at each component in the Vbe2 equation and see how it changes:
The change in Vbe1, delta_Vbe1 = delta_Ic/gm = delta_Ic * VT / Ic1
The change in Ic1*RL = delta_Ic * RL
To get delta_Vbe2 = 0 --> delta_IC * VT/ Ic1 = delta_IC * RL OR VT / Ic1 = RL
That's how you choose RL. Figure out your desired Ic1 & room temperature and that's it. You can also see how the voltage drop across RL ends up being VT like the linked page says. This would work perfectly if gm was constant over the full operating range, but alas, Ic1 changes and thus gm changes. To keep Ic2 constant, you'd need to change RL for every Ic1. So, in practice with a fixed RL, you get a "peaking" current source instead of a constant current source.
By choosing RL to satisfy delta_Vbe2 = 0, this will place you at the peak of the "peaking" graph.
Just picking out 2 random 2N3904's will mess up the accuracy. I don't really know by how much. The relationship of Ic1 and Ic2 at the peak being different by a factor of e can be calculated using the Ebers-Moll model. If you decide to try to figure this out, you can throw away the -1 factor in the Ebers-Moll equation since Vbe >> VT for this circuit's operating point.