;
; PURPOSE:    Calculate a field line from Rs to Rss using SS model
;             with Rss & Nmax (apar=0)
; WRITTEN:    05JAN2001	Xuepu Zhao 
; MODIFIED:   12MAR2003 Xuepu Zhao for adding rij
;
pro zoflij_ss,g,h,phri,thrj,np,xv,yv,zv,rv,tv,pv,$
              RIJ=rij,Rss=rss,NMAX=nmax,DSTEP=dstep,NSTEP=nstep
               

    if N_params() lt 1 then begin
       print,'pro iflij_ss,g,h,phri,thrj,np,xv,yv,zv,rv,tv,pv,$'
       print,'             rij=,rss=,apar=,nmax=,np='
       return
    endif

    if not keyword_set(rss) then rss=2.5
    if not keyword_set(nmax) then nmax=9
    if not keyword_set(dstep) then dstep=0.02
    if not keyword_set(nstep) then nstep=1500
    if not keyword_set(rij) then rij=1.0

    pi = !PI
    d2r=!DTOR
    xv=fltarr(nstep)
    yv=xv
    zv=xv
    rv=xv
    tv=xv
    pv=xv

    cth=COS(thrj)
    sth=SIN(thrj)
    cph=COS(phri)
    sph=SIN(phri)
    r=rij
    x = r*sth*cph
    y = r*sth*sph
    z = r*cth
    xv(0)=x
    yv(0)=y
    zv(0)=z
    rv(0)=rij
    tv(0)=thrj/!DTOR
    pv(0)=phri/!DTOR
    zpb_ss,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-1 do begin
        BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
        ds = sig * r * 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 1.0 OR r gt rss then begin
        ;  if sig eq 1 then np=(is-1) else np=-(is-1)
           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
           if z gt 0.0 then thr=atan(arg) else thr=atan(arg)+!pi
        endelse
        phr = atan(y,x)
        if phr lt 0 then phr=phr+2*!PI
        xv(is)=x
        yv(is)=y
        zv(is)=z
        rv(is)=r
        tv(is)=thr/d2r
        pv(is)=phr/d2r
        cth=COS(thr)
        sth=SIN(thr)
        cph=COS(phr)
        sph=SIN(phr)
        zpb_ss,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)
    tv=tv(0:is-1)
    pv=pv(0:is-1)
    return

end