***********************************************************************

       SUBROUTINE GUFLDT(XV,BTOS)

*  XV: 3 dim Position vector (input)
*  BTOS: 3 dim Field vector (output)  
*
*  created 12/15/98 R.Schicker 
*
***********************************************************************
	implicit none

#include "user.inc"
#include "fieldt.inc"

        real xv(3),btos(3)
        real xloc,yloc,zloc,rhog,rho2g,phig,phigd,phil,bz,br,bp
        real delz,delr,delp,delz1,delr1,delp1,delz2,delr2,delp2
        real dc,ds,wz,wr,wp,wzi,wri,wpi,w1,w2,w3,w4,w5,w6,w7,w8
        real dsec,dhsec,eps,radd,twopi,nsec,spol,dconv
        integer ifz,ifr,ifp 

        zloc=xv(3)
        if(zloc.lt.zflmin.or.zloc.gt.zflmax)go to 901

        xloc=xv(1)
        yloc=xv(2)

        rho2g=xloc**2+yloc**2
        rhog=sqrt(rho2g)

        if(rhog.gt.rflmax)go to 901

        eps=0.01
        radd=57.29577951
        twopi=6.283185308 
        dsec=60.
        dhsec=.5*dsec
        dconv=0.001

        if(rhog.gt.eps)then
           phig=atan2(yloc,xloc)
           if(phig.lt.0.)phig=phig+twopi
           phigd=radd*phig
           dc=cos(phig)
           ds=sin(phig)
         else
            rhog=0.
            phig=0.
            phigd=0.
            dc=0.
            ds=0.
         endif

         phil=mod(phigd,dsec)
         nsec=int((phigd-phil+eps)/dsec)+1
         if(phil.le.dhsec)spol=1.
         if(phil.gt.dhsec)then
            spol=-1.
            phil=dsec-phil
         endif

         delz=zloc-zflmin
         delr=rhog-rflmin
         delp=phil-pflmin

         delz2=mod(delz,zfldel)         
         delr2=mod(delr,rfldel)
         delp2=mod(delp,pfldel)

         delz1=delz-delz2
         delr1=delr-delr2
         delp1=delp-delp2

         wzi=delz2/zfldel
         wri=delr2/rfldel
         wpi=delp2/pfldel
         wz=1.-wzi
         wr=1.-wri
         wp=1.-wpi

         ifz=int((delz1+eps)/zfldel)+1
         ifr=int((delr1+eps)/rfldel)+1
         ifp=int((delp1+eps)/pfldel)+1

         w1=wz*wr*wp
         w2=wz*wri*wp
         w3=wz*wri*wpi
         w4=wz*wr*wpi
         w5=wzi*wr*wp
         w6=wzi*wri*wp
         w7=wzi*wri*wpi
         w8=wzi*wr*wpi

         bz=w1*tzfl(ifz,ifr,ifp)+w2*tzfl(ifz,ifr+1,ifp)
     1    +w3*tzfl(ifz,ifr+1,ifp+1)+w4*tzfl(ifz,ifr,ifp+1)   
     1    +w5*tzfl(ifz+1,ifr,ifp)+w6*tzfl(ifz+1,ifr+1,ifp)   
     1    +w7*tzfl(ifz+1,ifr+1,ifp+1)+w8*tzfl(ifz+1,ifr,ifp+1)

         br=w1*trfl(ifz,ifr,ifp)+w2*trfl(ifz,ifr+1,ifp)
     1    +w3*trfl(ifz,ifr+1,ifp+1)+w4*trfl(ifz,ifr,ifp+1)   
     1    +w5*trfl(ifz+1,ifr,ifp)+w6*trfl(ifz+1,ifr+1,ifp)   
     1    +w7*trfl(ifz+1,ifr+1,ifp+1)+w8*trfl(ifz+1,ifr,ifp+1)

         bp=w1*tpfl(ifz,ifr,ifp)+w2*tpfl(ifz,ifr+1,ifp)
     1    +w3*tpfl(ifz,ifr+1,ifp+1)+w4*tpfl(ifz,ifr,ifp+1)   
     1    +w5*tpfl(ifz+1,ifr,ifp)+w6*tpfl(ifz+1,ifr+1,ifp)   
     1    +w7*tpfl(ifz+1,ifr+1,ifp+1)+w8*tpfl(ifz+1,ifr,ifp+1)
         
         bz=-bz*fpol*spol*dconv
         br=-br*fpol*spol*dconv
         bp=-bp*fpol*dconv

         btos(1)=dc*br-ds*bp
         btos(2)=ds*br+dc*bp
         btos(3)=bz

         go to 999
 
 901     continue

         btos(1)=0.
         btos(2)=0.
         btos(3)=0.

 999     continue

       return
       end