add flux pumped JPA to README

README.md

@@ -53,7 +53,7 @@ Pkg.status()
 Generate a netlist using circuit components including capacitors `C`, inductors `L`, Josephson junctions described by the Josephson inductance `Lj`, mutual inductors described by the mutual coupling coefficient `K`, and resistors `R`. See the [SPICE netlist format](https://duckduckgo.com/?q=spice+netlist+format), docstrings, and examples below for usage. Run the harmonic balance analysis using `hbnlsolve` to solve a nonlinear system at one operating point, `hblinsolve` to solve a linear (or linearized) system at one or more frequencies, or `hbsolve` to run both analyses. Add a question mark `?` in front of a function to access the docstring.
 
 # Examples:
-## Josephson parametric amplifier
+## Josephson parametric amplifier (JPA)
 A driven nonlinear LC resonator.
 
 ```julia
@@ -95,6 +95,7 @@ plot(
             freqindex=:
         ),
     )),
+    label="JosephsonCircuits.jl",
     xlabel="Frequency (GHz)",
     ylabel="Gain (dB)",
 )
@@ -132,7 +133,7 @@ plot!(wswrspice/(2*pi*1e9),10*log10.(abs2.(S11)),
 ![JPA simulation with JosephsonCircuits.jl and WRspice](https://qce.mit.edu/JosephsonCircuits.jl/jpa_WRspice.png)
 
 
-## Double-pumped Josephson parametric amplifier
+## Double-pumped Josephson parametric amplifier (JPA)
 ```julia
 using JosephsonCircuits
 using Plots
@@ -205,6 +206,141 @@ plot!(wswrspice/(2*pi*1e9),10*log10.(abs2.(S11)),
 
 ![Double pumped JPA simulation with JosephsonCircuits.jl and WRspice](https://qce.mit.edu/JosephsonCircuits.jl/jpa_double_pumped_WRspice.png)
 
+## Flux-pumped Josephson parametric amplifier (JPA)
+Circuit and parameters from [here](https://doi.org/10.1063/1.2964182
+). Please note that three wave mixing (3WM) and flux-biasing are relatively untested, so you may encounter bugs. Please file issues or PRs.
+```julia
+using JosephsonCircuits
+using Plots
+
+@variables R Cc Cj Lj Cr Lr Ll Ldc K Lg
+circuit = [
+    ("P1","1","0",1),
+    ("R1","1","0",R),
+    # a very large inductor so the DC node flux of this node isn't floating
+    ("L0","1","0",Lg), 
+    ("C1","1","2",Cc),
+    ("L1","2","3",Lr),
+    ("C2","2","0",Cr),
+    ("Lj1","3","0",Lj),
+    ("Cj1","3","0",Cj),
+    ("L2","3","4",Ll),
+    ("Lj2","4","0",Lj),
+    ("Cj2","4","0",Cj),
+    ("L3","5","0",Ldc), 
+    ("K1","L2","L3",K),
+    # a port with a very large resistor so we can apply the bias across the port
+    ("P2","5","0",2),
+    ("R2","5","0",1000.0),
+] 
+
+circuitdefs = Dict(
+    Lj =>219.63e-12,
+    Lr =>0.4264e-9,
+    Lg =>100.0e-9,
+    Cc => 16.0e-15,
+    Cj => 10.0e-15, 
+    Cr => 0.4e-12,
+    R => 50.0, 
+    Ll => 34e-12, 
+    K => 0.999, # the inverse inductance matrix for K=1.0 diverges, so set K<1.0
+    Ldc => 0.74e-12,
+)
+
+ws = 2*pi*(9.7:0.0001:9.8)*1e9
+wp = (2*pi*19.50*1e9,)
+Ip = 0.7e-6
+Idc = 140.3e-6
+# add the DC bias and pump to port 2
+sourcespumpon = [(mode=(0,),port=2,current=Idc),(mode=(1,),port=2,current=Ip)]
+Npumpharmonics = (16,)
+Nmodulationharmonics = (8,)
+@time jpapumpon = hbsolve(ws, wp, sourcespumpon, Nmodulationharmonics,
+    Npumpharmonics, circuit, circuitdefs, dc = true, threewavemixing=true,fourwavemixing=true) # enable dc and three wave mixing
+
+
+plot(
+    jpapumpon.linearized.w/(2*pi*1e9),
+    10*log10.(abs2.(
+        jpapumpon.linearized.S(
+            outputmode=(0,),
+            outputport=1,
+            inputmode=(0,),
+            inputport=1,
+            freqindex=:
+        ),
+    )),
+    xlabel="Frequency (GHz)",
+    ylabel="Gain (dB)",
+    label="JosephsonCircuits.jl",
+)
+```
+
+```
+  0.015623 seconds (22.07 k allocations: 80.082 MiB)
+```
+
+and compare with WRspice
+```julia
+using XicTools_jll
+
+# simulate the JPA in WRSPICE
+wswrspice=2*pi*(9.7:0.005:9.8)*1e9
+n = JosephsonCircuits.exportnetlist(circuit,circuitdefs);
+input = JosephsonCircuits.wrspice_input_paramp(n.netlist,wswrspice,[0.0,wp[1]],[Idc,2*Ip],[(0,1)],[(0,5),(0,5)];trise=10e-9,tstop=600e-9);
+
+# @time output = JosephsonCircuits.spice_run(input,JosephsonCircuits.wrspice_cmd());
+@time output = JosephsonCircuits.spice_run(input,XicTools_jll.wrspice());
+S11,S21=JosephsonCircuits.wrspice_calcS_paramp(output,wswrspice,n.Nnodes);
+
+# plot the output
+plot!(wswrspice/(2*pi*1e9),10*log10.(abs2.(S11)),
+    label="WRspice",
+    seriestype=:scatter)
+```
+
+```
+283.557011 seconds (26.76 k allocations: 7.205 GiB, 0.66% gc time)
+```
+
+![Flux pumped JPA simulation with JosephsonCircuits.jl and WRspice](https://qce.mit.edu/JosephsonCircuits.jl/jpa_flux_pumped_WRspice.png)
+
+Simulate the JPA frequency as a function of DC bias current:
+```julia
+ws = 2*pi*(8.0:0.01:11.0)*1e9
+currentvals = (-20:0.1:20)*1e-5
+outvals = zeros(Complex{Float64},length(ws),length(currentvals))
+Ip=0.0
+
+Npumpharmonics = (1,)
+Nmodulationharmonics = (1,)
+
+@time for (k,Idc) in enumerate(currentvals)
+    sources = [
+          (mode=(0,),port=2,current=Idc),
+          (mode=(1,),port=2,current=Ip),
+      ]
+    sol = hbsolve(ws,wp,sources,Nmodulationharmonics, Npumpharmonics,
+        circuit, circuitdefs;dc=true,threewavemixing=true,fourwavemixing=true)
+    outvals[:,k]=sol.linearized.S((0,),1,(0,),1,:)
+end
+
+plot(
+    currentvals/(1e-3),
+    ws/(2*pi*1e9),
+    10*log10.(abs2.(outvals)),
+    seriestype=:heatmap,
+    xlabel="bias current (mA)",
+    ylabel="frequency (GHz)",
+    title="S11 (dB), pump off",
+)
+```
+
+```
+0.219279 seconds (3.27 M allocations: 639.981 MiB, 20.84% gc time)
+```
+
+![JPA frequency vs DC bias current](https://qce.mit.edu/JosephsonCircuits.jl/jpa_vs_bias_current.png)
 
 
 ## Josephson traveling wave parametric amplifier (JTWPA)