;
; Determin the strenth and SIGN of B between Rcp and Rss by downward integrating
;          from r to Rcp using HCCSSS model
;
pro pb_csss,g,h,gc,hc,apar,Nmax,Rcp,Rss,rr,thri,phrj,brrr,btrr,bprr

    d2r=!DTOR
    nao=!VALUES.F_NAN
    dstep=0.01
    nstep=500
; ***
    r=rr
    cth=cos(thri)
    sth=sin(thri)
    cph=cos(phrj)
    sph=sin(phrj)
    x = r*sth*cph
    y = r*sth*sph
    z = r*cth
    pb0_csss,gc,hc,apar,Nmax,Rcp,Rss,rr,cth,sth,phrj,Br,Bt,Bp
    brrr=Br & btrr=Bt & bprr=Bp
    BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
    sig = -1.
    r0=rr
    for is1=0,nstep do begin
        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 gt r0 then begin
;          print,'Br,..: ',brrr,btrr,bprr
;          print,'pb_csss: r0, r,thi,phj= ',r0,r,thri/d2r,phrj/d2r
           brrr=nao
           btrr=nao
           bprr=nao
           goto,lbl
        endif
        if abs(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
        if y gt 0.0 then phr=atan(y,x) else phr=2*!PI+atan(y,x)
        cth=COS(thr)
        sth=SIN(thr)
        cph=COS(phr)
        sph=SIN(phr)
        if r lt Rcp then begin
           pb_hc,g,h,r,cth,sth,phr,Br,Bt,Bp
           bsign=Br/ABS(Br)
           brrr=brrr*bsign
           btrr=btrr*bsign
           bprr=bprr*bsign
           goto, lbl
        endif
        r0=r
        Pb0_csss,gc,hc,apar,Nmax,Rcp,Rss,r,cth,sth,phr,Br,Bt,Bp
        BB = sqrt(Br*Br+Bt*Bt+Bp*Bp)
    endfor
 lbl:

end