;
; PURPOSE:    Calculate Flux Tube Expansion factor at phri,thrj using SS
; WRITTEN:    15NOV2002	  Xuepu Zhao 
; MODIFIED:   30APR2003	  Xuepu Zhao Adding keyword RBOTTOM & RSS
;
pro zfteij_ss,g,h,phri,thrj,fte,br1,phd1,thd1,NMAX=nmax,RSS=rss,RBOTTOM=rbottom

    if N_params() lt 1 then begin
       print,'fteij_ss,g,h,phri,thrj,fte,br1,phd1,thd1,nmax=,rss=,rbottom='
       print,'INPUT: g,h,phri,thrj'
       print,'Free parameters: nmax, rss, rbottom'
       print,'OUTPUT:fte,br1,phd1,thd1' 
       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(rbottom) then rbottom=1.0

    dstep=0.02
    nstep=600
    pi = !pi
    d2r=!DTOR
    nao=!values.f_nan

    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
    pb_ss,g,h,r,cth,sth,phri,Br,Bt,Bp,Nmax=nmax
    Brss = Br
    rar2=(rbottom/Rss)^2.0
    if Br gt 0.0 then sig = -1 else sig = 1
    for is=0,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 gt r0 then begin
           fte=nao
           phd1=nao
           thd1=nao
        ;  print,'fteij_ss stop at: r0,r,phd,thd=',r0,r,phd,thr/d2r
           return
        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
           cth=COS(thr)
           sth=SIN(thr)
           pb_ss,g,h,r,cth,sth,phr,Br,Bt,Bp,Nmax=nmax
           fte = (abs(Br)/Brss)*rar2
           br1=Br
           rbb=br1/Brss
           if rbb gt 0 then fte = (abs(br1)/Brss)*rar2 else fte=nao
           phd1=phd
           thd1=thr/d2r
           return
        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
        cth=COS(thr)
        sth=SIN(thr)
        cph=COS(phr)
        sph=SIN(phr)
        pb_ss,g,h,r,cth,sth,phr,Br,Bt,Bp,Nmax=nmax
    endfor

end