Added Lippman-Schwinger solver

doc/src/forces/.forces.exerinfo

@@ -113,7 +113,7 @@
             {'answer': '',
              'file': None,
              'hints': [],
-             'solution': '',
+             'solution': "The following Fortran 95 program solves the above Lippmann-Schwinger equation. Python and C++ codes will be added later. Discussions of the results will also be added.  \n!bc fcod\nC     *******************************************************\nC         Example program used to evaluate the \nC         T-matrix following Kowalski's method (eqs V88 & V89\nC         in Brown and Jackson)  \nC         for positive energies only\nC         The program is set up for S-waves only\nC         Coded by : Morten Hjorth-Jensen\nC         Language : FORTRAN 90\nC     *******************************************************\n\n\n\nC               ******************************\nC                 Def of global variables\nC               ******************************\n\n\n      MODULE constants\n         DOUBLE PRECISION , PUBLIC :: p_mass, hbarc\n         PARAMETER (p_mass =938.926D0, hbarc = 197.327D0)\n      END MODULE constants\n\n      MODULE mesh_variables           \n         INTEGER, PUBLIC :: n_rel\n         PARAMETER(n_rel=48)\n         DOUBLE PRECISION, ALLOCATABLE, PUBLIC :: ra(:), wra(:)\n      END MODULE  mesh_variables\n\n\nC               ******************************\nC                   Begin of main program\nC               ******************************\n\n\n      PROGRAM t_matrix\n      USE mesh_variables\n      IMPLICIT NONE\n      INTEGER istat\n\n      ALLOCATE( ra (n_rel), wra (n_rel),  STAT=istat )\n      CALL rel_mesh               ! rel mesh & weights\n      CALL t_channel              ! calculate the T-matrix\n      DEALLOCATE( ra,wra, STAT=istat )\n\n      END PROGRAM t_matrix\n\nC     *********************************************************\nC                    obtain the t-mtx\nC                    vkk is the box potential\nC                    f_mtx is equation V88 og Brown & Jackson\nC     *********************************************************\n\n      SUBROUTINE t_channel\n      USE mesh_variables\n      IMPLICIT NONE\n      INTEGER istat, i,j\n      DIMENSION vkk(:,:),f_mtx(:),t_mtx(:)\n      DOUBLE PRECISION, ALLOCATABLE :: vkk,t_mtx,f_mtx\n      DOUBLE PRECISION t_shell\n\n      ALLOCATE(vkk (n_rel,n_rel), STAT=istat)\n      CALL v_pot_yukawa(vkk) ! set up the box potential in routine vpot\n      ALLOCATE(t_mtx (n_rel), STAT=istat) ! allocate space in heap for T\n      ALLOCATE(f_mtx (n_rel), STAT=istat) ! allocate space for f\n      DO i=1,n_rel               ! loop over positive energies e=k^2\n         CALL f_mtx_eq(f_mtx,vkk,i)  ! solve eq. V88\n         CALL principal_value(vkk,f_mtx,i,t_shell) ! solve Eq. V89 \n         DO j=1,n_rel    ! the t-matrix\n            t_mtx(j)=f_mtx(j)*t_shell\n            IF(j == i) WRITE(6,*) ra(i) ,t_mtx(i)\nc     &                            DATAN(-ra(i)*t_mtx(i))\n         ENDDO\n      ENDDO\n      DEALLOCATE(vkk ,  STAT=istat)\n      DEALLOCATE(t_mtx, f_mtx,  STAT=istat)\n 1000 FORMAT( I3, 2F12.6) \n\n      END SUBROUTINE t_channel\n\nC     ***********************************************************\nC          The analytical expression for the box potential\nC          of exercise 1 and 12\nC          vkk is in units of fm^-2 (14 MeV/41.47Mevfm^2, where \nC          41.47= \\hbarc^2/mass_nucleon),  \nC          ra are mesh points in rel coordinates, units of fm^-1\nC     ***********************************************************\n\n      SUBROUTINE v_pot_box(vkk)\n      USE mesh_variables\n      USE constants\n      IMPLICIT NONE\n      INTEGER i,j\n      DOUBLE PRECISION  vkk, r_0, v_0, a, b, fac\n      PARAMETER(r_0=2.7d0,v_0=0.33759d0)  !r_0 in fm, v_0 in fm^-2 \n      DIMENSION vkk(n_rel,n_rel)\n\n      DO i=1,n_rel     ! set up the free potential\n         DO j=1,i-1 \n            a=ra(i)+ra(j)\n            b=ra(i)-ra(j)\n            fac=v_0/(2.d0*ra(i)*ra(j))\n            vkk(j,i)=fac*(DSIN(a*r_0)/a-DSIN(b*r_0)/b)\n            vkk(i,j)=vkk(j,i)\n         ENDDO\n         fac=v_0/(2.d0*(ra(i)**2))\n         vkk(i,i)=fac*(DSIN(2.d0*ra(i)*r_0)/(2.d0*ra(i))-r_0)\n      ENDDO\n\n      END  SUBROUTINE v_pot_box\n\nC     ***********************************************************\nC          The analytical expression for  a Yukawa potential\nC          in the l=0 channel\nC          vkk is in units of fm^-2,  \nC          ra are mesh points in rel coordinates, units of fm^-1\nC          The parameters here are those of the Reid-Soft core\nC          potential, see Brown and Jackson eq. A(4)\nC     ***********************************************************\n\n      SUBROUTINE v_pot_yukawa(vkk)\n      USE mesh_variables\n      USE constants\n      IMPLICIT NONE\n      INTEGER i,j\n      DOUBLE PRECISION  vkk, mu1, mu2, mu3, v_1, v_2, v_3, a, b, fac\n      PARAMETER(mu1=0.49d0,v_1=-0.252d0) \n      PARAMETER(mu2=7.84d0,v_2=-39.802d0) \n      PARAMETER(mu3=24.01d0,v_3=156.359d0) \n      DIMENSION vkk(n_rel,n_rel)\n\n      DO i=1,n_rel     ! set up the free potential\n         DO j=1,i \n            a=(ra(j)+ra(i))**2\n            b=(ra(j)-ra(i))**2\n            fac=1./(4.d0*ra(i)*ra(j))\n            vkk(j,i)=v_1*fac*DLOG((a+mu1)/(b+mu1))+\n     &               v_2*fac*DLOG((a+mu2)/(b+mu2))+\n     &               v_3*fac*DLOG((a+mu3)/(b+mu3))\n            vkk(i,j)=vkk(j,i)\n         ENDDO\n      ENDDO\n\n      END  SUBROUTINE v_pot_yukawa\n      \nC     **************************************************\nC         Solves eq. V88 \nC         and returns < p | f_mtx | n_pole =k>\nC     **************************************************\n\n      SUBROUTINE f_mtx_eq(f_mtx,vkk,n_pole)\n      USE mesh_variables\n      USE constants\n      IMPLICIT NONE\n      INTEGER i, j, int, istat, n_pole\n      DOUBLE PRECISION f_mtx,vkk,dp,deriv,pih,xsum\n      DIMENSION dp(1),deriv(1)\n      DIMENSION f_mtx(n_rel),vkk(n_rel,n_rel),a(:,:),fu(:),q(:),au(:)\n      DOUBLE PRECISION, ALLOCATABLE :: fu, q, au, a\n\n      pih=2.D0/ACOS(-1.D0)\n      ALLOCATE( a (n_rel,n_rel), STAT=istat)\n      DO i=1,n_rel\n         ALLOCATE(fu(n_rel), q(n_rel), au(n_rel), STAT=istat)\n         DO j=1,n_rel\n            fu(j)=vkk(i,j)-vkk(i,n_pole)*vkk(n_pole,j)/\n     &                 vkk(n_pole,n_pole)\n         ENDDO\n         DO j=1,n_rel\n            IF(j /= n_pole ) THEN     ! regular part\n               a(j,i)=pih*fu(j)*wra(j)*(ra(j)**2)/\n     &                (ra(j)**2-ra(n_pole)**2)\n            ELSEIF(j == n_pole) THEN  ! use l'Hopitals rule to get pole term\n               dp(1)=ra(j)             \n               CALL spls3(ra,fu,n_rel,dp,deriv(1),1,q,au,2,0) \n               a(j,i)=pih*wra(j)*ra(j)/2.d0*deriv(1)\n            ENDIF\n         ENDDO\n         DEALLOCATE(fu, q, au, STAT=istat)   ! free space in heap \n         a(i,i)=a(i,i)+1.D0\n      ENDDO\n      CALL matinv(a, n_rel)      ! Invert the matrix a\n      DO j=1,n_rel               ! multiply inverted matrix a with dim less pot\n         xsum=0.D0\n         DO i=1,n_rel\n            xsum=xsum+a(i,j)*vkk(i,n_pole)/vkk(n_pole,n_pole)  ! gives f-matrix in V88\n         ENDDO\n         f_mtx(j)=xsum\n      ENDDO\n      DEALLOCATE (a, STAT=istat)\n\n      END SUBROUTINE f_mtx_eq\n\nC     **************************************************\nC         Solves the principal value integral of V89\nC         returns the t-matrix for k=k, t_shell\nC     **************************************************\n\n      SUBROUTINE principal_value(vkk,f_mtx,n_pole,t_shell) \n      USE mesh_variables\n      IMPLICIT NONE\n      DOUBLE PRECISION vkk, f_mtx, t_shell, sum, pih, deriv, term\n      DIMENSION deriv(1)\n      DIMENSION vkk(n_rel, n_rel), f_mtx(n_rel),fu(:), q(:), au(:)\n      DOUBLE PRECISION, ALLOCATABLE :: fu, q, au\n      INTEGER n_pole, i, istat\n\n      ALLOCATE(fu(n_rel), q(n_rel), au(n_rel), STAT=istat)\n      sum=0.D0\n      pih=2.D0/ACOS(-1.D0)\n      DO i=1,n_rel\n         fu(i)=vkk(n_pole,i)*f_mtx(i)\n      ENDDO\n      DO i=1,n_pole-1  ! integrate up to the pole - 1 mesh\n         term=fu(i)*(ra(i)**2)-fu(n_pole)*(ra(n_pole)**2)\n         sum=sum+pih*wra(i)*term/(ra(i)**2-ra(n_pole)**2)\n      ENDDO       ! here comes the pole part\n      CALL spls3(ra,fu,n_rel,ra(n_pole),deriv,1,au,q,2,0)      \n      sum=sum+pih*wra(n_pole)*(fu(n_pole)+ra(n_pole)*deriv(1)/2.d0) \n      DO i=n_pole+1,n_rel  ! integrate from pole + 1mesh pt to infty\n         term=fu(i)*(ra(i)**2)-fu(n_pole)*(ra(n_pole)**2)\n         sum=sum+pih*wra(i)*term/(ra(i)**2-ra(n_pole)**2)\n      ENDDO\n      t_shell=vkk(n_pole,n_pole)/(1.d0+sum)\n      DEALLOCATE (fu, q, au, STAT=istat)\n\n      END SUBROUTINE principal_value\n\nC         ***********************************************\nC             Set up of relative mesh and weights\nC         ***********************************************\n\n      SUBROUTINE rel_mesh\n      USE mesh_variables\n      IMPLICIT NONE\n      INTEGER i\n      DOUBLE PRECISION pih,u,s,xx,c,h_max\n      PARAMETER (c=0.75)\n      DIMENSION u(n_rel), s(n_rel)\n\n      pih=ACOS(-1.D0)/2.D0\n      CALL gausslegendret (0.D0,1.d0,n_rel,u,s)\n      DO i=1,n_rel\n         xx=pih*u(i)\n         ra(i)=DTAN(xx)*c\n         wra(i)=pih*c/DCOS(xx)**2*s(i)\n      ENDDO\n\n      END SUBROUTINE rel_mesh\n\nC     *********************************************************\nC            Routines to do mtx inversion, from Numerical\nC            Recepies, Teukolsky et al. Routines included\nC            below are MATINV, LUDCMP and LUBKSB. See chap 2\nC            of Numerical Recepies for further details\nC            Recoded in FORTRAN 90 by M. Hjorth-Jensen\nC     *********************************************************\n\n      SUBROUTINE matinv(a,n)\n      IMPLICIT REAL*8(A-H,O-Z)\n      DIMENSION a(n,n)\n      INTEGER istat\n      DOUBLE PRECISION, ALLOCATABLE :: y(:,:)\n      INTEGER, ALLOCATABLE :: indx(:)\n\n      ALLOCATE (y( n, n), STAT =istat)\n      ALLOCATE ( indx (n), STAT =istat)\n      DO i=1,n\n         DO j=1,n\n            y(i,j)=0.\n         ENDDO\n         y(i,i)=1.\n      ENDDO\n      CALL  ludcmp(a,n,indx,d)\n      DO j=1,n\n         call lubksb(a,n,indx,y(1,j))\n      ENDDO\n      DO i=1,n\n         DO j=1,n\n            a(i,j)=y(i,j)\n         ENDDO\n      ENDDO\n      DEALLOCATE ( y, STAT=istat)\n      DEALLOCATE ( indx, STAT=istat)\n\n      END SUBROUTINE matinv\n \n\n      SUBROUTINE LUDCMP(A,N,INDX,D)\n      IMPLICIT REAL*8(A-H,O-Z)\n      PARAMETER (TINY=1.0E-20)\n      DIMENSION A(N,N),INDX(N)\n      INTEGER istat\n      DOUBLE PRECISION, ALLOCATABLE :: vv(:)\n\n      ALLOCATE ( vv(n), STAT = istat)\n      D=1.\n      DO I=1,N\n         AAMAX=0.\n         DO J=1,N\n            IF (ABS(A(I,J)) > AAMAX) AAMAX=ABS(A(I,J))\n         ENDDO\n         IF (AAMAX == 0.) PAUSE 'Singular matrix.'\n         VV(I)=1./AAMAX\n      ENDDO\n      DO J=1,N\n         IF (J > 1) THEN\n            DO I=1,J-1\n               SUM=A(I,J)\n               IF (I > 1)THEN\n                  DO K=1,I-1\n                     SUM=SUM-A(I,K)*A(K,J)\n                  ENDDO\n                  A(I,J)=SUM\n               ENDIF\n            ENDDO\n         ENDIF\n         AAMAX=0.\n         DO I=J,N\n            SUM=A(I,J)\n            IF (J > 1)THEN\n               DO K=1,J-1\n                  SUM=SUM-A(I,K)*A(K,J)\n               ENDDO\n               A(I,J)=SUM\n            ENDIF\n            DUM=VV(I)*ABS(SUM)\n            IF (DUM >= AAMAX) THEN\n               IMAX=I\n               AAMAX=DUM\n            ENDIF\n         ENDDO\n         IF (J /= IMAX)THEN\n            DO K=1,N\n               DUM=A(IMAX,K)\n               A(IMAX,K)=A(J,K)\n               A(J,K)=DUM\n            ENDDO\n            D=-D\n            VV(IMAX)=VV(J)\n         ENDIF\n         INDX(J)=IMAX\n         IF(J /= N)THEN\n            IF(A(J,J) == 0.) A(J,J)=TINY\n            DUM=1./A(J,J)\n            DO I=J+1,N\n               A(I,J)=A(I,J)*DUM\n            ENDDO\n         ENDIF\n      ENDDO\n      IF(A(N,N) == 0.)  A(N,N)=TINY\n      DEALLOCATE ( vv, STAT = istat)\n\n      END SUBROUTINE LUDCMP\n \n      SUBROUTINE LUBKSB(A,N,INDX,B)\n      implicit real*8(a-h,o-z)\n      DIMENSION A(N,N),INDX(N),B(N)\n      II=0\n      DO I=1,N\n         LL=INDX(I)\n         SUM=B(LL)\n         B(LL)=B(I)\n         IF (II /= 0)THEN\n            DO J=II,I-1\n               SUM=SUM-A(I,J)*B(J)\n            ENDDO\n         ELSE IF (SUM /= 0.) THEN\n            II=I\n         ENDIF\n         B(I)=SUM\n      ENDDO\n      DO I=N,1,-1\n         SUM=B(I)\n         IF (I < N)THEN\n            DO J=I+1,N\n               SUM=SUM-A(I,J)*B(J)\n            ENDDO\n         ENDIF\n         B(I)=SUM/A(I,I)\n      ENDDO\n\n      END SUBROUTINE lubksb\n\n!ec",
              'text': 'Write a small program which calculates the latter expression\nand use this potential to compute the $T$-matrix at positive energies\nfor $l=0$.\nCompare your results to those obtained with a  box potential given by\n!bt\n\\[\nV(r)=\\left\\{ \\begin{array}{cc} -V_0& r < R_0 \\\\\n                                0  & r > R_0 \\end{array} \\right.\n\\]\n\n!et\nMake a plot of the\ntwo $T$-matrices for energies up to 300 MeV in the lab frame\nand comment your results.\n\nFinally, a warning, the above central potential is fitted  \nto data from approximately \n20 MeV to some 300 MeV. This means that results outside\nthe data set should be taken seriously.'},
             {'answer': '',
              'file': None,

doc/src/forces/forces.aux

@@ -56,4 +56,5 @@
 \@writefile{toc}{\contentsline {paragraph}{b)}{31}{section*.21}}
 \@writefile{toc}{\contentsline {paragraph}{a)}{31}{section*.23}}
 \@writefile{toc}{\contentsline {paragraph}{b)}{31}{section*.24}}
-\@writefile{toc}{\contentsline {paragraph}{c)}{31}{section*.25}}
+\@writefile{toc}{\contentsline {paragraph}{Solution.}{31}{section*.25}}
+\@writefile{toc}{\contentsline {paragraph}{c)}{38}{section*.26}}

doc/src/forces/forces.do.txt

@@ -1811,6 +1811,393 @@ Finally, a warning, the above central potential is fitted
 to data from approximately 
 20 MeV to some 300 MeV. This means that results outside
 the data set should be taken seriously.
+!bsol
+The following Fortran 95 program solves the above Lippmann-Schwinger equation. Python and C++ codes will be added later. Discussions of the results will also be added.  
+!bc fcod
+C     *******************************************************
+C         Example program used to evaluate the 
+C         T-matrix following Kowalski's method (eqs V88 & V89
+C         in Brown and Jackson)  
+C         for positive energies only
+C         The program is set up for S-waves only
+C         Coded by : Morten Hjorth-Jensen
+C         Language : FORTRAN 90
+C     *******************************************************
+
+
+
+C               ******************************
+C                 Def of global variables
+C               ******************************
+
+
+      MODULE constants
+         DOUBLE PRECISION , PUBLIC :: p_mass, hbarc
+         PARAMETER (p_mass =938.926D0, hbarc = 197.327D0)
+      END MODULE constants
+
+      MODULE mesh_variables           
+         INTEGER, PUBLIC :: n_rel
+         PARAMETER(n_rel=48)
+         DOUBLE PRECISION, ALLOCATABLE, PUBLIC :: ra(:), wra(:)
+      END MODULE  mesh_variables
+
+
+C               ******************************
+C                   Begin of main program
+C               ******************************
+
+
+      PROGRAM t_matrix
+      USE mesh_variables
+      IMPLICIT NONE
+      INTEGER istat
+
+      ALLOCATE( ra (n_rel), wra (n_rel),  STAT=istat )
+      CALL rel_mesh               ! rel mesh & weights
+      CALL t_channel              ! calculate the T-matrix
+      DEALLOCATE( ra,wra, STAT=istat )
+
+      END PROGRAM t_matrix
+
+C     *********************************************************
+C                    obtain the t-mtx
+C                    vkk is the box potential
+C                    f_mtx is equation V88 og Brown & Jackson
+C     *********************************************************
+
+      SUBROUTINE t_channel
+      USE mesh_variables
+      IMPLICIT NONE
+      INTEGER istat, i,j
+      DIMENSION vkk(:,:),f_mtx(:),t_mtx(:)
+      DOUBLE PRECISION, ALLOCATABLE :: vkk,t_mtx,f_mtx
+      DOUBLE PRECISION t_shell
+
+      ALLOCATE(vkk (n_rel,n_rel), STAT=istat)
+      CALL v_pot_yukawa(vkk) ! set up the box potential in routine vpot
+      ALLOCATE(t_mtx (n_rel), STAT=istat) ! allocate space in heap for T
+      ALLOCATE(f_mtx (n_rel), STAT=istat) ! allocate space for f
+      DO i=1,n_rel               ! loop over positive energies e=k^2
+         CALL f_mtx_eq(f_mtx,vkk,i)  ! solve eq. V88
+         CALL principal_value(vkk,f_mtx,i,t_shell) ! solve Eq. V89 
+         DO j=1,n_rel    ! the t-matrix
+            t_mtx(j)=f_mtx(j)*t_shell
+            IF(j == i) WRITE(6,*) ra(i) ,t_mtx(i)
+c     &                            DATAN(-ra(i)*t_mtx(i))
+         ENDDO
+      ENDDO
+      DEALLOCATE(vkk ,  STAT=istat)
+      DEALLOCATE(t_mtx, f_mtx,  STAT=istat)
+ 1000 FORMAT( I3, 2F12.6) 
+
+      END SUBROUTINE t_channel
+
+C     ***********************************************************
+C          The analytical expression for the box potential
+C          of exercise 1 and 12
+C          vkk is in units of fm^-2 (14 MeV/41.47Mevfm^2, where 
+C          41.47= \hbarc^2/mass_nucleon),  
+C          ra are mesh points in rel coordinates, units of fm^-1
+C     ***********************************************************
+
+      SUBROUTINE v_pot_box(vkk)
+      USE mesh_variables
+      USE constants
+      IMPLICIT NONE
+      INTEGER i,j
+      DOUBLE PRECISION  vkk, r_0, v_0, a, b, fac
+      PARAMETER(r_0=2.7d0,v_0=0.33759d0)  !r_0 in fm, v_0 in fm^-2 
+      DIMENSION vkk(n_rel,n_rel)
+
+      DO i=1,n_rel     ! set up the free potential
+         DO j=1,i-1 
+            a=ra(i)+ra(j)
+            b=ra(i)-ra(j)
+            fac=v_0/(2.d0*ra(i)*ra(j))
+            vkk(j,i)=fac*(DSIN(a*r_0)/a-DSIN(b*r_0)/b)
+            vkk(i,j)=vkk(j,i)
+         ENDDO
+         fac=v_0/(2.d0*(ra(i)**2))
+         vkk(i,i)=fac*(DSIN(2.d0*ra(i)*r_0)/(2.d0*ra(i))-r_0)
+      ENDDO
+
+      END  SUBROUTINE v_pot_box
+
+C     ***********************************************************
+C          The analytical expression for  a Yukawa potential
+C          in the l=0 channel
+C          vkk is in units of fm^-2,  
+C          ra are mesh points in rel coordinates, units of fm^-1
+C          The parameters here are those of the Reid-Soft core
+C          potential, see Brown and Jackson eq. A(4)
+C     ***********************************************************
+
+      SUBROUTINE v_pot_yukawa(vkk)
+      USE mesh_variables
+      USE constants
+      IMPLICIT NONE
+      INTEGER i,j
+      DOUBLE PRECISION  vkk, mu1, mu2, mu3, v_1, v_2, v_3, a, b, fac
+      PARAMETER(mu1=0.49d0,v_1=-0.252d0) 
+      PARAMETER(mu2=7.84d0,v_2=-39.802d0) 
+      PARAMETER(mu3=24.01d0,v_3=156.359d0) 
+      DIMENSION vkk(n_rel,n_rel)
+
+      DO i=1,n_rel     ! set up the free potential
+         DO j=1,i 
+            a=(ra(j)+ra(i))**2
+            b=(ra(j)-ra(i))**2
+            fac=1./(4.d0*ra(i)*ra(j))
+            vkk(j,i)=v_1*fac*DLOG((a+mu1)/(b+mu1))+
+     &               v_2*fac*DLOG((a+mu2)/(b+mu2))+
+     &               v_3*fac*DLOG((a+mu3)/(b+mu3))
+            vkk(i,j)=vkk(j,i)
+         ENDDO
+      ENDDO
+
+      END  SUBROUTINE v_pot_yukawa
+
+C     **************************************************
+C         Solves eq. V88 
+C         and returns < p | f_mtx | n_pole =k>
+C     **************************************************
+
+      SUBROUTINE f_mtx_eq(f_mtx,vkk,n_pole)
+      USE mesh_variables
+      USE constants
+      IMPLICIT NONE
+      INTEGER i, j, int, istat, n_pole
+      DOUBLE PRECISION f_mtx,vkk,dp,deriv,pih,xsum
+      DIMENSION dp(1),deriv(1)
+      DIMENSION f_mtx(n_rel),vkk(n_rel,n_rel),a(:,:),fu(:),q(:),au(:)
+      DOUBLE PRECISION, ALLOCATABLE :: fu, q, au, a
+
+      pih=2.D0/ACOS(-1.D0)
+      ALLOCATE( a (n_rel,n_rel), STAT=istat)
+      DO i=1,n_rel
+         ALLOCATE(fu(n_rel), q(n_rel), au(n_rel), STAT=istat)
+         DO j=1,n_rel
+            fu(j)=vkk(i,j)-vkk(i,n_pole)*vkk(n_pole,j)/
+     &                 vkk(n_pole,n_pole)
+         ENDDO
+         DO j=1,n_rel
+            IF(j /= n_pole ) THEN     ! regular part
+               a(j,i)=pih*fu(j)*wra(j)*(ra(j)**2)/
+     &                (ra(j)**2-ra(n_pole)**2)
+            ELSEIF(j == n_pole) THEN  ! use l'Hopitals rule to get pole term
+               dp(1)=ra(j)             
+               CALL spls3(ra,fu,n_rel,dp,deriv(1),1,q,au,2,0) 
+               a(j,i)=pih*wra(j)*ra(j)/2.d0*deriv(1)
+            ENDIF
+         ENDDO
+         DEALLOCATE(fu, q, au, STAT=istat)   ! free space in heap 
+         a(i,i)=a(i,i)+1.D0
+      ENDDO
+      CALL matinv(a, n_rel)      ! Invert the matrix a
+      DO j=1,n_rel               ! multiply inverted matrix a with dim less pot
+         xsum=0.D0
+         DO i=1,n_rel
+            xsum=xsum+a(i,j)*vkk(i,n_pole)/vkk(n_pole,n_pole)  ! gives f-matrix in V88
+         ENDDO
+         f_mtx(j)=xsum
+      ENDDO
+      DEALLOCATE (a, STAT=istat)
+
+      END SUBROUTINE f_mtx_eq
+
+C     **************************************************
+C         Solves the principal value integral of V89
+C         returns the t-matrix for k=k, t_shell
+C     **************************************************
+
+      SUBROUTINE principal_value(vkk,f_mtx,n_pole,t_shell) 
+      USE mesh_variables
+      IMPLICIT NONE
+      DOUBLE PRECISION vkk, f_mtx, t_shell, sum, pih, deriv, term
+      DIMENSION deriv(1)
+      DIMENSION vkk(n_rel, n_rel), f_mtx(n_rel),fu(:), q(:), au(:)
+      DOUBLE PRECISION, ALLOCATABLE :: fu, q, au
+      INTEGER n_pole, i, istat
+
+      ALLOCATE(fu(n_rel), q(n_rel), au(n_rel), STAT=istat)
+      sum=0.D0
+      pih=2.D0/ACOS(-1.D0)
+      DO i=1,n_rel
+         fu(i)=vkk(n_pole,i)*f_mtx(i)
+      ENDDO
+      DO i=1,n_pole-1  ! integrate up to the pole - 1 mesh
+         term=fu(i)*(ra(i)**2)-fu(n_pole)*(ra(n_pole)**2)
+         sum=sum+pih*wra(i)*term/(ra(i)**2-ra(n_pole)**2)
+      ENDDO       ! here comes the pole part
+      CALL spls3(ra,fu,n_rel,ra(n_pole),deriv,1,au,q,2,0)      
+      sum=sum+pih*wra(n_pole)*(fu(n_pole)+ra(n_pole)*deriv(1)/2.d0) 
+      DO i=n_pole+1,n_rel  ! integrate from pole + 1mesh pt to infty
+         term=fu(i)*(ra(i)**2)-fu(n_pole)*(ra(n_pole)**2)
+         sum=sum+pih*wra(i)*term/(ra(i)**2-ra(n_pole)**2)
+      ENDDO
+      t_shell=vkk(n_pole,n_pole)/(1.d0+sum)
+      DEALLOCATE (fu, q, au, STAT=istat)
+
+      END SUBROUTINE principal_value
+
+C         ***********************************************
+C             Set up of relative mesh and weights
+C         ***********************************************
+
+      SUBROUTINE rel_mesh
+      USE mesh_variables
+      IMPLICIT NONE
+      INTEGER i
+      DOUBLE PRECISION pih,u,s,xx,c,h_max
+      PARAMETER (c=0.75)
+      DIMENSION u(n_rel), s(n_rel)
+
+      pih=ACOS(-1.D0)/2.D0
+      CALL gausslegendret (0.D0,1.d0,n_rel,u,s)
+      DO i=1,n_rel
+         xx=pih*u(i)
+         ra(i)=DTAN(xx)*c
+         wra(i)=pih*c/DCOS(xx)**2*s(i)
+      ENDDO
+
+      END SUBROUTINE rel_mesh
+
+C     *********************************************************
+C            Routines to do mtx inversion, from Numerical
+C            Recepies, Teukolsky et al. Routines included
+C            below are MATINV, LUDCMP and LUBKSB. See chap 2
+C            of Numerical Recepies for further details
+C            Recoded in FORTRAN 90 by M. Hjorth-Jensen
+C     *********************************************************
+
+      SUBROUTINE matinv(a,n)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION a(n,n)
+      INTEGER istat
+      DOUBLE PRECISION, ALLOCATABLE :: y(:,:)
+      INTEGER, ALLOCATABLE :: indx(:)
+
+      ALLOCATE (y( n, n), STAT =istat)
+      ALLOCATE ( indx (n), STAT =istat)
+      DO i=1,n
+         DO j=1,n
+            y(i,j)=0.
+         ENDDO
+         y(i,i)=1.
+      ENDDO
+      CALL  ludcmp(a,n,indx,d)
+      DO j=1,n
+         call lubksb(a,n,indx,y(1,j))
+      ENDDO
+      DO i=1,n
+         DO j=1,n
+            a(i,j)=y(i,j)
+         ENDDO
+      ENDDO
+      DEALLOCATE ( y, STAT=istat)
+      DEALLOCATE ( indx, STAT=istat)
+
+      END SUBROUTINE matinv
+
+
+      SUBROUTINE LUDCMP(A,N,INDX,D)
+      IMPLICIT REAL*8(A-H,O-Z)
+      PARAMETER (TINY=1.0E-20)
+      DIMENSION A(N,N),INDX(N)
+      INTEGER istat
+      DOUBLE PRECISION, ALLOCATABLE :: vv(:)
+
+      ALLOCATE ( vv(n), STAT = istat)
+      D=1.
+      DO I=1,N
+         AAMAX=0.
+         DO J=1,N
+            IF (ABS(A(I,J)) > AAMAX) AAMAX=ABS(A(I,J))
+         ENDDO
+         IF (AAMAX == 0.) PAUSE 'Singular matrix.'
+         VV(I)=1./AAMAX
+      ENDDO
+      DO J=1,N
+         IF (J > 1) THEN
+            DO I=1,J-1
+               SUM=A(I,J)
+               IF (I > 1)THEN
+                  DO K=1,I-1
+                     SUM=SUM-A(I,K)*A(K,J)
+                  ENDDO
+                  A(I,J)=SUM
+               ENDIF
+            ENDDO
+         ENDIF
+         AAMAX=0.
+         DO I=J,N
+            SUM=A(I,J)
+            IF (J > 1)THEN
+               DO K=1,J-1
+                  SUM=SUM-A(I,K)*A(K,J)
+               ENDDO
+               A(I,J)=SUM
+            ENDIF
+            DUM=VV(I)*ABS(SUM)
+            IF (DUM >= AAMAX) THEN
+               IMAX=I
+               AAMAX=DUM
+            ENDIF
+         ENDDO
+         IF (J /= IMAX)THEN
+            DO K=1,N
+               DUM=A(IMAX,K)
+               A(IMAX,K)=A(J,K)
+               A(J,K)=DUM
+            ENDDO
+            D=-D
+            VV(IMAX)=VV(J)
+         ENDIF
+         INDX(J)=IMAX
+         IF(J /= N)THEN
+            IF(A(J,J) == 0.) A(J,J)=TINY
+            DUM=1./A(J,J)
+            DO I=J+1,N
+               A(I,J)=A(I,J)*DUM
+            ENDDO
+         ENDIF
+      ENDDO
+      IF(A(N,N) == 0.)  A(N,N)=TINY
+      DEALLOCATE ( vv, STAT = istat)
+
+      END SUBROUTINE LUDCMP
+
+      SUBROUTINE LUBKSB(A,N,INDX,B)
+      implicit real*8(a-h,o-z)
+      DIMENSION A(N,N),INDX(N),B(N)
+      II=0
+      DO I=1,N
+         LL=INDX(I)
+         SUM=B(LL)
+         B(LL)=B(I)
+         IF (II /= 0)THEN
+            DO J=II,I-1
+               SUM=SUM-A(I,J)*B(J)
+            ENDDO
+         ELSE IF (SUM /= 0.) THEN
+            II=I
+         ENDIF
+         B(I)=SUM
+      ENDDO
+      DO I=N,1,-1
+         SUM=B(I)
+         IF (I < N)THEN
+            DO J=I+1,N
+               SUM=SUM-A(I,J)*B(J)
+            ENDDO
+         ENDIF
+         B(I)=SUM/A(I,I)
+      ENDDO
+
+      END SUBROUTINE lubksb
+!ec
+!esol
 !esubex
 !bsubex
 The parameters of the box potential are chosen to fit a