Example-Josephson-junction.md

-All the files (inputs and outputs) for this example are accessible via [Zenodo](https://zenodo.org/records/18404682).
+All input and output files for this example are available on [Zenodo](https://zenodo.org/records/18404682).
+
+[[_TOC_]]
 
 # Target
-Simulate quasi-1D Josephson dynamics. Initialize the simulation of a unitary Fermi gas confined in an external potential
+Simulate quasi-1D Josephson dynamics in a unitary Fermi gas confined by
 ```math
 V^{(ext)}_{\sigma}(x) = \frac{1}{2}\omega_x^2 x^2 + V_0e^{-2x^2/w^2} + k x,
 ```
-where the first term corresponds to harmonic confinement, the second describes the potential barrier separating the cloud into two reservoirs, and the last introduces a small linear tilt to generate a population imbalance between the left and right reservoirs.
+* Harmonic trap: confinement.
+* Gaussian barrier: splits the cloud into two reservoirs.
+* Linear tilt: generates an initial population imbalance.
 
-# Step 1: Creating working folders
-We start by copying template folders for static (st) and time-dependent (td) simulations:
+# Step 1: Create Working Directories
+Copy the project templates:
 ```bash
 cp -r $WSLDA/st-project-template ./st-jj-example
 cp -r $WSLDA/td-project-template ./td-jj-example
 ```
 
-# Step 2: Editing files for the static calculations
-All edits should be done in the working folder `st-jj-example`.
+# Step 2: Static Calculation Setup (`st-jj-example`)
 
 ## 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:
+Define lattice and functional:
 ```c
 /**
  * Define lattice size and lattice spacing.
@@ -34,8 +37,10 @@ Select the lattice size and functional, and inspect all other tags that can be r
 #define SLDA_FORCE_A1
 // ...
 ```
+Adjust other tags if needed.
+
 ## 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. 
+Define the external potential:
 ```c
 #include "../extensions/wslda_utils.h"
 // ...
@@ -50,13 +55,16 @@ double v_ext(int ix, int iy, int iz, int it, int spin, double *params, size_t ex
 }
 // ...
 ```
-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
+### Human-Readable Parameters
+
+Instead of specifying $`\omega_x`$  directly, define the Thomas–Fermi radius:
 ```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. 
+Barrier parameters are given in units of $`\varepsilon_F`$ and $`k_F`$.  
+Conversion is handled in `process_params(...)`:
 ```c
 void process_params(double *params, double *kF, double *mu, size_t extra_data_size, void *extra_data)
 {
@@ -72,14 +80,13 @@ void process_params(double *params, double *kF, double *mu, size_t extra_data_si
     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 logger form. 
+In this example, we use the default logger form. 
 
 ## 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 
+Key user parameters:
 ```bash
 # -------------- USER DEFINED PARAMETERS --------------
 # Data flow: [Read params from input file] -> [execute process_params( )] 
@@ -92,7 +99,7 @@ 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
+Other input parameters relevant for this run:
 ```bash
 # ...
 # ------------------- INITIALIZATION ------------------
@@ -125,7 +132,7 @@ init0maxiter  1000   # maximum number of iterations when searching for a  unifor
 # ...
 ``` 
 
-# Step 3: Compiling and running the static code
+# Step 3: Compile & Run (Static)
 To generate the binary for a quasi-1D simulation, execute the command:
 ```bash
 make 1d
@@ -143,8 +150,8 @@ mpirun -np 4 ./st-wslda-1d input.txt
 ```
 See [here](./example-jj/st/jj-r1.stdout) for full stdout. The file produced by the logger is [here](./example-jj/st/jj-r1.wlog).
 
-# Step 4: Analyzing the results of the static run
-As a simple example, we will plot the density profile and compute the population imbalance, defined as $`z=(N_L-N_R)/(N_L+N_R)`$, where $`N_L`$ and $`N_R`$ are particle numbers in left and right reservoirs, using Python 
+# Step 4: Analyze Static Results
+Example: density profile and population imbalance $`z=(N_L-N_R)/(N_L+N_R)`$, where $`N_L`$ and $`N_R`$ are particle numbers in left and right reservoirs
 ```python
 import numpy as np
 import matplotlib.pyplot as plt
@@ -169,9 +176,8 @@ The result is
 ![st-example-result](uploads/854d1a0926ca45a5c9891a618623d012/st-example-result.png)
 For a more advanced script for plotting static results, see [Zenodo](https://zenodo.org/records/18404682).
 
-# Step 5: Editing files for the time-dependent calculations
-All edits should be done in the working folder `td-jj-example`.
-The workflow is analogous to that for static calculations. Many of the functions can be copied & pasted from static files. 
+# Step 5: Time-Dependent Setup (`td-jj-example`)
+Workflow is analogous.
 
 ## Edit [predefines.h](./example-jj/td/predefines.h)
 Set compilation flags to be compatible with your problem and with data produced by static calcs. 
@@ -180,7 +186,7 @@ Set compilation flags to be compatible with your problem and with data produced
 Edit the function that defines the external potential. You can also add time dependence to your external potential function if you wish. 
 
 ## Edit [logger.h](./example-jj/td/logger.h)
-In this file user can customize the output. In this example, we use the default logger form. 
+In this example, we use the default logger form. 
 
 ## Edit [input.txt](./example-jj/td/input.txt)
 In the input file, the key parameters for this example are:
@@ -210,8 +216,8 @@ selfstart      1 # use Taylor expansion for first 5 steps?
 # ...
 ```
 
-# Step 6: Compiling and running the time-dependent code
-We do this part in an analogous way to the case of static calculations. Note that the time-dependent run requires GPUs to run. 
+# Step 6: Compile & Run (Time-Dependent)
+**GPU resources required.**  
 
 To generate the binary for a quasi-1D simulation, execute the command:
 ```bash
@@ -240,8 +246,33 @@ mpirun -np 1 ./td-wslda-1d input.txt
 ```
 See [here](./example-jj/td/jj-r2.stdout) for full stdout. The file produced by the logger is [here](./example-jj/td/jj-r2.wlog).
 
-# Step 7: Analyzing the results of the time-dependent run
-As an example of analysis, we will plot population imbalance and phase differnce accors the junction as a function of time, using simple python script
+# Step 7: Analyze Time Evolution
+Plot population imbalance and phase difference as a function of time.
 ```python
-# todo
+import numpy as np
+import matplotlib.pyplot as plt
+# to get wdata lib use: pip install wdata
+from wdata.io import WData, Var
+
+# Load data
+data = WData.load("./td-jj-example/jj-r2.wtxt")
+x=data.xyz[0]
+rho=data.rho_a+data.rho_b # total density for all time frames..
+delta=data.delta/data.eF  # normalized order parameter for all time frames..
+
+# compute initial population imbalance for each time frame...
+NL = np.sum(rho[:,x<0], axis=1)
+NR = np.sum(rho[:,x>0], axis=1)
+z=(NL-NR)/(NL+NR) # population imbalance
+# ... and phase difference
+phase_diff = np.array([np.angle(delta[frame, 128 + 10] / delta[frame, 128 - 10]) / np.pi for frame in range(delta.shape[0])])
+
+# plot population imbalance and phase difference as a function of time
+# make two panels
+fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6))
+ax1.plot(z, color='blue')
+ax2.plot(phase_diff, color='red')
+plt.show()
 ```
+![td-example-result](uploads/81c30bb26def96ce1487d38c950e8100/td-example-result.png)
+For a more advanced script for plotting time-dependent results, see [Zenodo](https://zenodo.org/records/18404682).
\ No newline at end of file