anchor.diff:                write(out,my_fmt,advance="no"), &
                 tranfreq,abscoef,jrot(ilevelf),energyf, jrot(ileveli),energyi 
                 !
  -              write(out,"(a4)",advance="no"), " <- "
  +              !write(out,"(a4)",advance="no"), " <- "
                 !
                 do kitem = 1,maxitems 
                   !
  
positive.code:              !
               do ifilter = 1,Nfilters
                 !
                 if (upper(ifilter)%mask/=trim(quantum_numbers(upper(ifilter)%i,ilevelf)).or.lower(ifilter)%mask/=trim(quantum_numbers(lower(ifilter)%i,ileveli))) cycle loop_tran
                 !
               enddo
               !
  
negative.code:        !
         return
      end subroutine gauher
  
  
      subroutine gauher_half(x,w,n)
         !
         ! modified for w = exp(-x^2/2)/sqrt(2 pi)
         !
         integer(ik) :: n,MAXIT
         real(rk) :: w(n),x(n)
         real(rk) ::  EPS,PIM4
         parameter (EPS=3.D-16,PIM4=.75112554446494248285870300477622_rk,MAXIT=40)
         integer(ik) :: i,its,j,m
         real(rk) ::p1,p2,p3,pp,z,z1
  
         !
         !if (n==20) then 
         !  !
         !     x(1:20)  = &
         !       (/-7.619048541679757_rk,-6.510590157013656_rk,-5.578738805893203_rk,&
         !       -4.734581334046057_rk,-3.943967350657318_rk,-3.18901481655339_rk,&
         !       -2.458663611172367_rk,-1.745247320814127_rk,-1.042945348802751_rk,&
         !       -0.346964157081356_rk, 0.346964157081356_rk, 1.042945348802751_rk,&
         !        1.745247320814127_rk, 2.458663611172367_rk, 3.18901481655339_rk,&
         !        3.943967350657316_rk, 4.734581334046057_rk, 5.578738805893202_rk,&
         !        6.510590157013653_rk, 7.619048541679757_rk/)
         !   w(1:20) =  (/0.000000000000126_rk, 0.000000000248206_rk, 0.000000061274903_rk,&
         !        0.00000440212109_rk, 0.000128826279962_rk, 0.00183010313108_rk,&
         !        0.013997837447101_rk, 0.061506372063977_rk, 0.161739333984_rk,&
         !        0.260793063449555_rk, 0.260793063449555_rk, 0.161739333984_rk,&
         !        0.061506372063977_rk, 0.013997837447101_rk, 0.00183010313108_rk,&
         !        0.000128826279962_rk, 0.00000440212109_rk, 0.000000061274903_rk,&
         !        0.000000000248206_rk, 0.000000000000126_rk/)
         !   return
         !  !
         !endif 
         !
         ! do only half 
         ! m=(n+1)/2
         !
         m=n
         !
         do i=1,m
           !
           if(i.eq.1)then
             z=sqrt(float(2*n+1))-1.85575*(2*n+1)**(-.16667)
           else if(i.eq.2)then
             z=z-1.14*n**.426/z
           else if (i.eq.3)then
             z=1.86*z-.86*x(1)
           else if (i.eq.4)then
             z=1.91*z-.91*x(2)
           else
             z=2.*z-x(i-2)
           endif
           do its=1,MAXIT
             p1=PIM4
             p2=0.0_rk
             do j=1,n
               p3=p2
               p2=p1
               p1=z*sqrt(2.0_rk/j)*p2-sqrt(real(j-1,ark)/real(j,ark))*p3
             enddo
             pp=sqrt(2.0_rk*n)*p2
             z1=z
             z=z1-p1/pp
             if(abs(z-z1).le.EPS) exit
             !
           enddo
           if (its==MAXIT) stop 'too many iterations in gauher'
           x(i)=z
           !
           ! exclude the symmetric point
           !x(n+1-i)=-z
           !
           w(i)=2.0_rk/(pp*pp)
           !
           ! exclude the symmetric point
           !w(n+1-i)=w(i)
           !
         enddo
         !
         w = w / sqrt(pi) * 2.0_rk
         x = x*sqrt(2.0_rk)
         !
         return
      end subroutine gauher_half
  
  
  
      !C  (C) Copr. 1986-92 Numerical Recipes Software 
  
    subroutine gauher_ark(xr,wr,n)
  