diff --git a/src/dispersion.f90 b/src/dispersion.f90
index 8ef76ea..1403a97 100644
--- a/src/dispersion.f90
+++ b/src/dispersion.f90
@@ -10,32 +10,43 @@ module dispersion
   real(wp_), dimension(npts+1), save :: ttv,extv
 !
 contains
-!
-!
-subroutine colddisp(xg,yg,npl,nprf,sox)  
-! solution cold dispersion relation
-!
+
+
+pure function colddisp(xg, yg, npl, sox) result(npr)
+  ! Cold plasma dispersion relation
+
+  ! Given the parallel component of refractive index returns
+  ! the orthogonal one.
+  !
+  ! Reference: IFP-CNR Internal Report FP 05/1 - App. A
+  !
   implicit none
-! arguments
-!  xg   = omegap^2/omega^2
-!  yg   = omegac/omega
-!  npl  = N parallel to B
-!  nprf = N perpendicular to B (cold)
-!  sox  = mode (-1 => O mode, +1 => X mode)
-  real(wp_) :: xg,yg,npl,nprf,sox
-! local variables
-  real(wp_) :: yg2,npl2,dnl,del,dxg
-!  
-  npl2=npl*npl
-  dnl=one-npl2
-  dxg=one-xg
-  yg2=yg*yg
-  del=sqrt(dnl*dnl+4.0_wp_*npl2*dxg/yg2)
-  nprf=sqrt(dxg-npl2-xg*yg2*(one+npl2+sox*del)/(dxg-yg2)/2.0_wp_)
-!
-end subroutine colddisp
-!
-!
+
+  ! subroutine arguments
+
+  ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
+  real(wp_), intent(in) :: xg, yg
+
+  ! cold refractive index: N∥, N⊥
+  real(wp_), intent(in)  :: npl
+  real(wp_)              :: npr
+
+  ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
+  real(wp_), intent(in) :: sox
+
+  ! local variables
+  real(wp_) :: yg2, npl2, dnl, del, dxg
+
+  npl2 = npl**2
+  yg2  = yg**2
+  dnl  = one - npl2
+  dxg  = one - xg
+  del  = sqrt(dnl**2 + 4.0_wp_*npl2*dxg/yg2)
+  npr  = sqrt(dxg - npl2 - xg*yg2 * (one + npl2 + sox*del)/(dxg - yg2)/2.0_wp_)
+
+end function colddisp
+
+
 !
 subroutine harmnumber(yg,mu,npl,nhmin,nhmax,iwr)
 ! computation of minimum and maximum harmonic
@@ -114,26 +125,35 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
   complex(wp_) :: cc0,cc2,cc4,discr,rdiscr,npra2,npra,npr,npr2,e330,e11,e22, &
                   e12,e13,e23,a11,a22,a33,a12,a13,a23,a330,a1122,a123,a330n, &
                   cc0t,cc2t,cc4t
-  complex(wp_) :: epsl(3,3,lrm),sepsl(3,3)
+
+  ! eij: dielectric tensor
+
+  complex(wp_) :: epsl(3,3,lrm)  !
+  complex(wp_) :: sepsl(3,3)     !  
+
 !
   err=0
   errnpr=one
-  npra2=nprf**2
+  npra2=nprf**2  ! first guest N⊥ (usually cold)
   npl2=npl*npl
   dnl=one-npl2
   imxx=abs(imx)
-!
-  if (fast.eq.1) then
+
+  ! Compute dielectric tensor elements
+  ! returns ε^l/X
+  if (fast == 1) then
     call diel_tens_wr(yg,mu,npl,a330,epsl,lrm)
   else
     call diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast)
   end if
-  a330= xg*a330
-  e330= one + a330
-!
-  do
-    do i=1,imxx
-!  
+
+  a330 = xg*a330
+  e330 = one + a330
+
+  tries_loop: do
+    iterations_loop: do i=1, imxx
+
+      ! sum the ε^l up to lrm
       do j=1,3
         do k=1,3
           sepsl(k,j)=czero
@@ -142,86 +162,84 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
           end do
         end do
       end do
-!
+
       npra=sqrt(npra2)
-!
+
       a11 = xg*sepsl(1,1)
       a22 = xg*sepsl(2,2)
       a12 = xg*sepsl(1,2)
       a33 = xg*sepsl(3,3)
       a13 = xg*sepsl(1,3)
       a23 = xg*sepsl(2,3)
-!      a31 = a13
-!      a32 =-a23
-!
+      ! a31 = a13
+      ! a32 =-a23
+
       e11 = one + a11
       e22 = one + a22
       e12 = a12
-!      e33 = e330 + npra2*a33
+      ! e33 = e330 + npra2*a33
       e13 = npra*a13
       e23 = npra*a23
-!      e21 =-e12
-!      e31 = e13
-!      e32 =-e23
-!
-!      if(i.gt.2.and.errnpr.lt.1.0e-3_wp_) exit
-!
-      cc4 = (a11 + dnl)*(one - a33) + (a13 + npl)*(a13 + npl)
-      cc2 =-a12*a12*(one - a33) + two*a23*a12*(a13 + npl) &
+      ! e21 =-e12
+      ! e31 = e13
+      ! e32 =-e23
+
+      ! A, B, C coefficients of the biquadratic (eq. 29)
+      ! AN⊥⁴ + BN⊥² + C = 0
+      cc4 = (a11 + dnl)*(one - a33) + (a13 + npl)*(a13 + npl)   ! A
+      cc2 =-a12*a12*(one - a33) + two*a23*a12*(a13 + npl) &     ! B
            -(e330 + (a22 + dnl)*(one - a33) - a23*a23)*(a11 + dnl) &
            -(a13 + npl)*(a13 + npl)*(a22 + dnl)
-      cc0 = e330*((a11 + dnl)*(a22 + dnl) + a12*a12)
-!
+      cc0 = e330*((a11 + dnl)*(a22 + dnl) + a12*a12)            ! C
+
       cc4t = cc4 - one
       cc2t = half*cc2 + dnl
       cc0t = cc0 - dnl*dnl
-!
+
       a1122 = a11*a22 + a12*a12
       a330n = a330 + dnl*a33
       a123  = a12*a23 - a13*a22
-!
+
+      ! Discriminant Δ = B² - 4*A*C
+      ! Note: this has been rewritten to avoid floating points cancellation
       discr = (cc2t*cc2t - cc0t*cc4t) &
              -( npl2*a1122 + dnl*a22*a330n + (dnl*a23)**2 + two*dnl*npl*a123 &
                + a1122*a330n + dnl*a13*a123 + dnl*a23*(a12*a13 + a11*a23) )
-!
-!      if(yg.gt.one) then
-!        s=sox
-!        if(dimag(discr).LE.zero) s=-s
-!      else
-!        s=-sox
-!        if(dimag(discr).ge.zero) s=-s
-!      end if
-!
+
       rdiscr=sqrt(discr)
-      if(dimag(rdiscr/cc4).gt.0.0d0) then
-!      if(dimag(discr).gt.0.0d0) then
-        s=sox
+
+      ! Choice of the solution (±) for the given mode
+      if(dimag(rdiscr/cc4) > zero) then
+        s = sox
       else
-        s=-sox
+        s = -sox
       end if
-!
-      npr2=(s*rdiscr - half*cc2)/cc4
-!
+
+      npr2 = (s*rdiscr - half*cc2)/cc4
+
       errnpr=abs(one-abs(npr2)/abs(npra2))
       if(i.gt.1.and.errnpr.lt.1.0e-5_wp_) exit
       npra2=npr2
-    end do
+    end do iterations_loop
 
     if(i.gt.imxx.and.imxx.gt.1) then
       if (imx.lt.0) then
+        ! first try failed
         err=1
         imxx=1
         npra2=nprf**2
       else
+        ! first try failed, no retry
         err=2
         exit
       end if
     else
+      ! converged or second try
       exit
     end if
-        
-  end do 
-!
+
+  end do tries_loop
+
   if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).ge.huge(one).or. &
      abs(npr2).le.tiny(one)) then
     npr2=czero