anchor.diff:    do 160 l=1,il
       if(abs(el(l)).le.pieni) goto 160
       kz1=kz(l)+1
  -!       goto(20,40,80,60,40,60,100,100,140),kz1
  -!       goto 160
  -!Eric
   !-----------------------------------------------------------------------
   !  DRIFTLENGTH
   !-----------------------------------------------------------------------
  -    if (kz1.eq.1) then
  +    if(kz1 == 1) then
         goto 20
   !-----------------------------------------------------------------------
   !  RECTANGULAR MAGNET
   !  HORIZONTAL
   !-----------------------------------------------------------------------
  -    elseif (kz1.eq.2.or.kz1.eq.5) then
  -40     fokm=el(l)*ed(l)
  -      if(abs(fokm).le.pieni) goto 20
  -      if(kz1.eq.2) then
  +    elseif(kz1 == 2 .or. kz1 == 5) then
  +40    fokm=el(l)*ed(l)
  +      if(abs(fokm) <= pieni) goto 20
  +      if(kz1 == 2) then
           ih1=1
           ih2=2
         else
  
positive.code:         as(4,ih2,j,l) = ((-one*two)*as6)*fok1
          as(5,ih2,j,l) = (as6*fok1)*fok1
          as(6,ih2,j,l) = as6
  50     continue
        goto 160
      elseif (kz1.eq.4.or.kz1.eq.6) then
  !-----------------------------------------------------------------------
  !  SEKTORMAGNET
  !  HORIZONTAL
  !-----------------------------------------------------------------------
  60     fokm=el(l)*ed(l)
        if(abs(fokm).le.pieni) goto 20
        if(kz1.eq.4) then
          ih1=1
          ih2=2
        else
  !  SECTOR MAGNET VERTICAL
          ih1=2
          ih2=1
        endif
        do 70 j=1,napx
          fok  = fokm/dpsq(j)
          rho  = (one/ed(l))*dpsq(j)
          si   = sin_mb(fok)
          co   = cos_mb(fok)
          rhoc = (rho*(one-co))/dpsq(j)
          siq=si/dpsq(j)
          al(1,ih1,j,l)=co
          al(2,ih1,j,l)=rho*si
          al(3,ih1,j,l)=(-one*si)/rho                            !hr01
          al(4,ih1,j,l)=co
          al(5,ih1,j,l)=((-one*dpsv(j))*rhoc)*c1e3                  !hr01
          al(6,ih1,j,l)=((-one*dpsv(j))*siq)*c1e3                   !hr01
          sm12=el(l)-al(1,ih1,j,l)*al(2,ih1,j,l)
          sm23=al(2,ih1,j,l)*al(3,ih1,j,l)
          as(1,ih1,j,l)=(el(l)*(one-rvv(j))-rvv(j)*((dpsv(j)**2/(four*dpd(j)))*sm12+dpsv(j)*(el(l)-al(2,ih1,j,l))))*c1e3
          as(2,ih1,j,l)=(-one*rvv(j))*((dpsv(j)/(two*rho*dpsq(j)))*sm12-dpd(j)*siq)
          as(3,ih1,j,l)=(-one*rvv(j))*(((dpsv(j)*rho)/(two*dpsq(j)))*sm23-dpd(j)*rhoc)
          as(4,ih1,j,l)=((-one*rvv(j))*sm23)/c2e3                   !hr01
          as(5,ih1,j,l)=((-one*rvv(j))*sm12)/((c4e3*rho)*rho) !hr01
          as(6,ih1,j,l)=((-one*rvv(j))*(el(l)+al(1,ih1,j,l)*al(2,ih1,j,l)))/c4e3
  !--VERTIKAL
          as(6,ih2,j,l)=((-one*rvv(j))*al(2,ih2,j,l))/c2e3             !hr01
  70     continue
        goto 160
      elseif (kz1.eq.3) then
  !-----------------------------------------------------------------------
  !  QUADRUPOLE
  !  FOCUSSING
  !-----------------------------------------------------------------------
  80   do 90 j=1,napx
          fok=ekv(j,l)*oidpsv(j)
          aek(j)=abs(fok)
          hi(j)=sqrt(aek(j))
          fi(j)=el(l)*hi(j)
  
negative.code:         ih1=2
          ih2=1
        endif
        do j=1,napx
          fok  = fokm/dpsq(j)
          rho  = (one/ed(l))*dpsq(j)
          fok1 = (tan_mb(fok*half))/rho
  