;
; PURPOSE:    Calculate a field line from Rss to Rbottom, expansion factor or
;             foot points using HCSS model with Rss & Nmax
; WRITTEN:    08MAY2003	Xuepu Zhao 
; MODIFIED:   fte=(BB/Brss)/(Rss/Rbottom)^2
;
pro ziflij_ss,g,h,phri,thrj,np,xv,yv,zv,rv,tv,pv,fte,RSS=rss,NMAX=nmax,$
              DSTEP=dstep,NSTEP=nstep,RBOTTOM=rbottom

    if N_params() lt 1 then begin
       print,'pro ziflij_ss,g,h,phri,thrj,np,xv,yv,zv,rv,tv,pv,fte,'
       print,'              rss=,nmax=,dstep=,nstep=,rbottom='
       return
    endif

    if not keyword_set(nmax) then begin
       sz=SIZE(g)
       nmax=sz(1)-1
    endif
    if not keyword_set(rss) then rss=2.5
    if not keyword_set(dstep) then dstep=0.02
    if not keyword_set(nstep) then nstep=600
    if not keyword_set(rbottom) then rbottom=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=Rss
    x = r*sth*cph
    y = r*sth*sph
    z = r*cth
    xv(0)=x
    yv(0)=y
    zv(0)=z
    rv(0)=Rss
    tv(0)=thrj/d2r
    pv(0)=phri/d2r
    zpb_ss,g,h,r,cth,sth,phri,Br,Bt,Bp,Nmax=nmax
    Brss = Br
    rar=rbottom/rss
    rar2=rar*rar
    if Br gt 0.0 then sig = -1 else sig = 1
    for is=1,nstep do begin
        BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
        r0=r
        x0=x & y0=y & z0=z
        ds = sig * r0 * 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 = x0 + dx
        y = y0 + dy
        z = z0 + dz
        r = sqrt(x*x + y*y + z*z)
        if r ge r0 then begin
           fte=0
           np=0
           goto, lbl
        endif
        if r lt Rbottom then begin
           dstep_end=r0-rbottom
           ds = sig * r0 * dstep_end /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 = x0 + dx
           y = y0 + dy
           z = z0 + dz
           r = rbottom
           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 gt 0.0 then phd=phr/d2r else phd=phr/d2r+360.0
           xv(is)=x
           yv(is)=y
           zv(is)=z
           rv(is)=r
           tv(is)=thr/d2r
           pv(is)=phd
           cth=COS(thr)
           sth=SIN(thr)
           zpb_ss,g,h,r,cth,sth,phr,Br,Bt,Bp,Nmax=nmax
           fte = (abs(Br)/Brss)*rar2
           np=is+1
           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 gt 0.0 then phd=phr/d2r else phd=phr/d2r+360.0
        xv(is)=x
        yv(is)=y
        zv(is)=z
        rv(is)=r
        tv(is)=thr/d2r
        pv(is)=phd
        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
 lbl: return

end