;
; PURPOSE:    Calculate a field line from Rs to Rss using SS model
;             with Rss & Nmax (apar=0)
; WRITTEN:    05JAN2001 Xuepu Zhao
; MODIFIED:   26JAN2006 for tstep, ...
;
pro oflijss_gaha06,g,h,phri,thrj,np,xv,yv,zv,rv,tvv,pv,$
               Rss=rss,RBOTTOM=rbottom,DSTEP=dstep

    if N_params() lt 1 then begin
       print,'pro oflijss_gaha06,g,h,phri,thrj,np,xv,yv,zv,rv,tvv,pv,$'
       print,'               rss=,dstep=,nstep='
       return
    endif

    if not keyword_set(rss) then rss=2.5
    if not keyword_set(rbottom) then rbottom=1.0
    if not keyword_set(dstep) then dstep=0.01
    nstep=1000
    sz=SIZE(g) & nmax=sz(1)-1
    pi = !PI
    d2r=!DTOR
    xv=fltarr(nstep)
    yv=xv
    zv=xv
    rv=xv
    tvv=xv
    pv=xv

    rv(0)=1.0
    tvv(0)=thrj/!DTOR
    pv(0)=phri/!DTOR
    cph=COS(phri)
    sph=SIN(phri)
    cth=COS(thrj)
    sth=SIN(thrj)
    r=Rbottom
    x = r*sth*cph
    y = r*sth*sph
    z = r*cth
    xv(0)=x
    yv(0)=y
    zv(0)=z
    pb_ss06,g,h,r,cth,sth,phri,Br,Bt,Bp,nmax=nmax
    if Br gt 0.0 then sig = 1 else sig = -1
    for is=1,nstep do begin           ; not (nstep-1) for finding something
        BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
        ds = sig * r * dstep / BB
    ;   ds = sig * dstep / BB
        dx = ds*(Br*sth*cph+Bt*cth*cph-Bp*sph)
        dy = ds*(Br*sth*sph+Bt*cth*sph+Bp*cph)
        dz = ds*(Br*cth-Bt*sth)
        x = x + dx
        y = y + dy
        z = z + dz
        r = sqrt(x*x + y*y + z*z)
        if r lt Rbottom OR r gt rss then begin
           if sig eq 1 then np=is else np=-is
           goto, lbl
        endif
        if z eq 0.0 then begin
           thr=!pi/2.0
        endif else begin
           arg=sqrt(x*x +y*y)/z
           thr=atan(arg)
           if z lt 0.0 then thr=thr+!pi
        endelse
        phr = atan(y,x)
        if phr lt 0.0 then phr=phr+2*!PI
        xv(is)=x
        yv(is)=y
        zv(is)=z
        rv(is)=r
        tvv(is)=thr/d2r
        pv(is)=phr/d2r
        cth=COS(thr)
        sth=SIN(thr)
        cph=COS(phr)
        sph=SIN(phr)
        pb_ss06,g,h,r,cth,sth,phr,Br,Bt,Bp,nmax=nmax
    endfor
    np=0
 lbl:
    xv=xv(0:is-1)
    yv=yv(0:is-1)
    zv=zv(0:is-1)
    rv=rv(0:is-1)
    tvv=tvv(0:is-1)
    pv=pv(0:is-1)
    return

end