;
; PURPOSE:    Calculate a field line from Rss to Rbottom, expansion factor and/or
;             foot points using HCSS model with Rss & Nmax
; WRITTEN:    14JUN2000 Xuepu Zhao
; MODIFIED:   16JUN2005 Xuepu Zhao for moving (r lt Rbottom) after pb_ss
; MODIFIED:   01FEB2006 for tstep, ...
;
pro iflijss_gaha06,g,h,phri,thrj,np,fte,xv,yv,zv,rv,tvv,pv,Rss=rss,$
               RBOTTOM=Rbottom,DSTEP=dstep

    if N_params() lt 1 then begin
       print,'pro iflijss_gaha06,g,h,phri,thrj,np,fte,xv,yv,zv,rv,tvv,pv,'
       print,'               rss=,rbottom=,dstep='        
       return
    endif
;
    if not keyword_set(rss) then rss=2.5           ; Added 19APR2006
    if not keyword_set(Rbottom) then Rbottom=1.00
    if not keyword_set(dstep) then dstep=0.02
    nstep=500
    sz=SIZE(g)
    nmax=sz(1)-1
    d2r=!DTOR
    xv=fltarr(nstep)
    yv=xv
    zv=xv
    rv=xv
    tvv=xv
    pv=xv
;
    rv(0)=rss
    tvv(0)=thrj/d2r
    pv(0)=phri/d2r
    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
;
    pb_ss06,g,h,r,cth,sth,phri,Br,Bt,Bp,Nmax=nmax
    BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
    brss = Br
    rar=rss/Rbottom
    rar2=rar*rar
    r0=r
    if Br gt 0.0 then sig = -1 else sig = 1
    for is=1,nstep do begin
    ;   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 ge r0 then begin
           fte=0
           np=0
           goto, lbl
        endif
        if (r-Rbottom) lt 0 then begin
           fte = (ABS(Br)/brss)/rar2
           np=is*brss/ABS(brss)
           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
        tvv(is)=thr/d2r
        pv(is)=phd
        cth=COS(thr)
        sth=SIN(thr)
        pb_ss06,g,h,r,cth,sth,phr,Br,Bt,Bp,Nmax=nmax
        BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
        r0=r
    endfor
 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