Example-Josephson-junction.md

+All the files (inputs and outputs) for this example are accessible via [Zenodo](https://zenodo.org/records/18404682).
+
 # Target
 Simulate quasi-1D Josephson dynamics. Initialize the simulation of a unitary Fermi gas confined in an external potential
 ```math
@@ -13,7 +15,7 @@ cp -r $WSLDA/td-project-template ./td-jj-example
 ```
 
 # Step 2: Editing files for the static calculations
-## Edit predefines.h
+## Edit [predefines.h](./example-jj/st/predefines.h)
 Select the lattice size and functional, and inspect all other tags that can be relevant for your problem:
 ```c
 /**
@@ -25,8 +27,116 @@ Select the lattice size and functional, and inspect all other tags that can be r
 // ...
 #define FUNCTIONAL SLDAE
 /**
- * Sets the effective mass to be equal 1.
+ * Sets the effective mass to be equal to 1.
  * */
 #define SLDA_FORCE_A1
 // ...
-```
+```
\ No newline at end of file
+## Edit [problem-definition.h](./example-jj/st/problem-definition.h)
+In this file, you provide a formula for the external potential. We will parametrize the function so it can be controlled via an input file. 
+```c
+#include "../extensions/wslda_utils.h"
+// ...
+double v_ext(int ix, int iy, int iz, int it, int spin, double *params, size_t extra_data_size, void *extra_data)
+{
+    double x = DX*(ix-NX/2);
+
+    double V_ho = harmonic_oscillator_smooth_edges(x, params[11], params[2]*LX*0.5, params[3]*LX*0.5, 1.0);
+    double V_barrier = params[4] * exp(-2.0*pow(x/params[5],2));
+    double V_tilt = params[6]*x/(LX/2); // linear tilt potential
+    return V_ho + V_barrier + V_tilt;
+}
+// ...
+```
+It is more convenient to provide parameter values that define the potential in human-readable units. For example, rather than specifying the trapping frequency $`\omega_x`$ directly, user provides the desired Thomas-Fermi radius $`x_{TF}`$, which is related to the chemical potential through
+```math
+    \frac{1}{2}\omega_x^2 x_{TF}^2 = \mu
+    \quad \Rightarrow \quad
+    \omega_x = \frac{\sqrt{2\mu}}{x_{TF}}.
+```
+Another example concerns the barrier height $`V_0`$ and width $`w`$, which are more conveniently expressed in terms of the Fermi energy $`\varepsilon_F`$ and Fermi wave vector $`k_F`$. To do a conversion from human-readable units to code units, we can use the `process_params(...)` function. 
+```c
+void process_params(double *params, double *kF, double *mu, size_t extra_data_size, void *extra_data)
+{
+    // if akF is active, then overwrite the value of input->sclgth
+    // value of kF is generated by the function referencekF(...)
+    if(input->akF!=0.0) input->sclgth = input->akF/kF[0];
+
+    // PROCESS INPUT FILE PARAMETERS
+    // Definition of R_TF:  0.5*m*omega^2 * R_TF^2 = mu
+    //                      => omega = sqrt(2*mu/m)/R_TF
+    params[11] = sqrt(2.*mu[SPINA])/ (params[1]*LX*0.5); // omega_x
+    double eF = kF[0]*kF[0]/2.0;
+    params[4] = params[4]*eF; // convert barrier height to internal units
+    params[5] = params[5]/kF[0]; // convert barrier width to internal units
+}
+
+```
+
+## Edit [logger.h](./example-jj/st/logger.h)
+In this file user can customize the output. In this example, we use the default form of logger. 
+
+## Edit [input.txt](./example-jj/st/input.txt)
+The input file provides run-time parameters. The most important in this case are the parameters that define the potential 
+```bash
+# -------------- USER DEFINED PARAMETERS --------------
+# Data flow: [Read params from input file] -> [execute process_params( )] 
+#                    -> [pass params to functions]
+params[1] = 0.80  # Thomas-Fermi radius for x direction, in units of LX/2
+params[2] = 0.90  # Smoothing parameter x1 for harmonic_oscillator_smooth_edges
+params[3] = 0.97  # Smoothing parameter x2 for harmonic_oscillator_smooth_edges
+params[4] = 0.24  # height of the barrier, in eF units
+params[5] = 5.0   # width of barrier, in 1/kF units
+params[6] = 0.015 # tilt potential (value at the edge of the trap)
+# ...
+```
+Other input parameters relevant for this run are
+```bash
+# ...
+# ------------------- INITIALIZATION ------------------
+inittype     0       # create uniform solution and start from it,
+# ...
+# ------------------- INPUT/OUTPUT ------------------
+outprefix    jj-r1   # all output files will start with this prefix
+writewf      1       # write wf at the end of computation: yes=1, no=0
+                     # we will evolve them by td-wslda-1d code
+# variables to write
+writevar    rho delta j V_ext 
+sclgth      -10000.0 # scattering length in units of lattice spacing
+referencekF 1.0      # in this case, I can define the value of kF=1
+# ...
+# ----------------- SELF-CONSISTENT LOOP -----------------
+energyconveps 1.0e-6 # energy convergence epsilon
+npartconveps  1.0e+6 # large number means that this condition will always be satisfied
+muchange      0.0    # do not change chemical potential, 
+                     # so computation will be for fixed chemical potential
+maxiters     200     # maximum number of iterations
+spinsymmetry 1       # to impose Na=Nb
+# ...
+# ----------- IN CASE OF: inittype=0 ----------------
+# See: Wiki -> Initialization of the solver
+init0na      0.01689 # required density for component a
+init0nb      0.01689 # and component b, so corresponding kF=1
+init0muchange 0.5    # change rate of chemical potentials, default muchange
+init0eps      1.0e-9 # epsilon for convergence for uniform solution
+init0maxiter  1000   # maximum number of iterations when searching for a  uniform solution
+# ...
+``` 
+
+# Step 3: Compiling and running the static code
+To generate the binary for a quasi-1D simulation, execute the command:
+```bash
+make 1d
+```
+If the compilation completes successfully, a binary named `st-wslda-1d` will be created.
+This program should be executed within an MPI environment.
+```bash
+mpirun -np 4 ./st-wslda-1d input.txt
+# START OF THE MAIN FUNCTION
+# CODE: ST-WSLDA-1D
+...
+# ALGORITHM CONVERGED!
+...
+# EXTRA SAVING ITERATION DONE.
+```
+See [here](./example-jj/st/jj-r1.stdout) for full stdout.