diff --git a/CMakeLists.txt b/CMakeLists.txt
index da9391c..0a8f2e4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,6 +15,7 @@ set(SOURCES
   src/main.cpp
   src/mandelbrot.s
   src/mandelbrotshader.cpp
+  src/julia.s
   src/juliashader.cpp
   src/utilities.cpp
   # ...
diff --git a/README.md b/README.md
index eb1c02c..d4a9652 100644
--- a/README.md
+++ b/README.md
@@ -1,22 +1,24 @@
-# A simple shader for AZUR aiming at computating the Mandelbrot or Julia set
+# A simple shader for AZUR aiming at computing the Mandelbrot and the Julia fractal sets
 
-This small demo is based on Azur.
-It uses ultra fast rendering pipeline with specific shaders :
- - Mandelbrot shader written in ASM by Lephe
- - Julia shader written in C with the libnum fast fixed point arithmetic lib
+## What is this addin doing ?  
 
+This small demo is based on [Azur](https://gitea.planet-casio.com/Lephenixnoir/Azur). It uses its ultra fast rendering pipeline with specific shaders :
+- Mandelbrot shader written in ASM by Lephe'
+- Julia shader written in C with the libnum fast fixed point arithmetic lib
 
-Keys are :
- - [EXIT] to leave to the Operating System
- - [F1] to [F4] to show/switch off some informations on the screen :
-   o [F1] : hide everything
-   o [F2] : minimal output FPS + rendering time
-   o [F3] : detailed update/rendering time
-   o [F4] : memory usage
- - [F5]/[F6] : switch fractal mode
-   o [F5] : set Mandelbrot fractal
-   o [F6] : set Julia fractal
+Both Mandelbrot - *named as per Benoît MANDELBROT (1924-2010)* - and Julia - *named as per Gaston JULIA (1893-1978)* - fractals are based on the convergence determination of the equation `Z_(n+1) = Z_n^2 + C`.
 
+##  What can we do ?
 
-Have fun !!!
+Control keys are :
+[**EXIT**] to leave to the Operating System
+[**F1**] to [**F4**] to show/switch off some informations on the screen :
+ - [ ] [F1] : hide everything
+ - [ ] [F2] : minimal output FPS + rendering time
+ - [ ] [F3] : detailed update/rendering time
+ - [ ] [F4] : memory usage
+[**F5**]/[**F6**] : switch fractal mode
+ - [ ] [F5] : set Mandelbrot fractal
+ - [ ] [F6] : set Julia fractal
 
+  Have fun !!!
diff --git a/src/julia.s b/src/julia.s
new file mode 100644
index 0000000..68ececd
--- /dev/null
+++ b/src/julia.s
@@ -0,0 +1,76 @@
+
+/*
+int julia(int32_t cr, int32_t ci, uint32_t t_squared, int steps, int32_t zr, int32_t zi);
+(Pour moi : t_squared = (2 * 2) << 16, steps = 50)
+Si renvoie 0 -> dans la fractale
+Sinon renvoie un nombre entre 1...steps indiquant la "distance" (pour le dégradé)
+*/
+
+.global _julia
+.text
+
+/* Iteration of z = z²+c with threshold checks
+   @r4   Real part of c, cr
+   @r5   Imaginary part of c, ci
+   @r6   Squared threshold t²
+   @r7   Number of steps
+   If c is in the set, returns 0. Otherwise, returns the remaining number of
+   iterations to be checked before the process finised. */
+_julia:
+	mov.l	r8, @-r15
+	mov.l	r9, @-r15
+	mov.l	r10, @-r15
+	mov.l	r11, @-r15
+
+	/* Start with z=z²=0 (z=x+iy) */
+	mov	r8, r0
+	mov	r9, r1
+	mov	#0, r2
+	mov	#0, r3
+
+.loop:
+	/* Start xy */
+	dmuls.l	r0, r1
+
+	/* Update x = x²-y²+cr */
+	mov	r2, r0
+	sub	r3, r0
+	add	r4, r0
+
+	/* Get y = 2xy+ci while starting x²
+	   TODO: Hope there's no overflow in that [shll r1]. It might be
+	   TODO: provable because 2xy overflowing probably implies a lower
+	   TODO: bound on x²+y² at the previous iteration. */
+	sts	mach, r10
+	sts	macl, r1
+	dmuls.l	r0, r0
+	xtrct	r10, r1
+	shll	r1
+	add	r5, r1
+
+	/* Get x² and keep the integer part at full precision for ∥z∥²; at the
+	   same time, start y² */
+	sts	mach, r10
+	sts	macl, r2
+	dmuls.l	r1, r1
+	xtrct	r10, r2
+
+	/* Get y² and compute ∥z∥² = x²+y² in 32:0 format */
+	sts	mach, r11
+	sts	macl, r3
+	add	r11, r10
+	xtrct	r11, r3
+
+	/* Compare ∥z∥² to t², return if the threshold is reached */
+	cmp/ge	r6, r10
+	bt	.end
+
+	/* Continue looping */
+	dt	r7
+	bf	.loop
+
+.end:
+	mov.l	@r15+, r11
+	mov.l	@r15+, r10
+	rts
+	mov	r7, r0
\ No newline at end of file
diff --git a/src/juliashader.cpp b/src/juliashader.cpp
index a09f2b1..10be8fb 100644
--- a/src/juliashader.cpp
+++ b/src/juliashader.cpp
@@ -83,6 +83,9 @@ void azrp_julia( void )
     prof_leave(azrp_perf_cmdgen);
 }
 
+extern "C" {
+    extern int julia(int32_t cr, int32_t ci, uint32_t t_squared, int steps, int32_t zr, int32_t zi);
+}
 
 void azrp_shader_julia( void *uniforms, void *comnd, void *fragment )
 {
@@ -100,18 +103,21 @@ void azrp_shader_julia( void *uniforms, void *comnd, void *fragment )
     libnum::num zrsave = cmd->xmin;
     libnum::num zisave = cmd->ymax - libnum::num(cmd->current_frag*16)*cmd->yinc;
 
+    int u;
+
     int offset = 0;
 
     for( int j=0; j<16; j++)
     {      
         for( int i=0; i<azrp_width; i++ )
         {
+
             zr = zrsave;
             zi = zisave;
             zr2 = zr * zr;
             zi2 = zi * zi;
 
-            int u=0;
+            u=0;
 
             while((zr2+zi2)<4 && (u<maxiter))
             {
@@ -167,7 +173,6 @@ void azrp_shader_julia( void *uniforms, void *comnd, void *fragment )
         zisave -= cmd->yinc;
 
         offset += azrp_width;
-
     }
 
     cmd->current_frag++;