;
; Determin the strenth and SIGN of B at rr between Rcp and Rss
; by downward integrating from rr to Rcp using HCCSSS model
; NOTE: brrr, ... in mT
; MODIFIED: 08DEC2005	Xuepu Zhao
;
pro ppb_csss,ga,ha,gc,hc,rr,thri,phrj,APAR=apar,NMAX=nmax,RCP=rcp,RSS=rss,brrr,btrr,bprr

    if N_params() lt 1 then begin
       print,'pb_csss,ga,ha,gc,hc,rr,thri,phrj,apar=,Nmax=,Rcp=,Rss=,'
       print,'        brrr,btrr,bprr'
       print,'OUTPUT: brrr,btrr,bprr'
       return
    endif
    if not keyword_set(apar) then apar=0.2
    if not keyword_set(nmax) then nmax=15
    if not keyword_set(rcp) then rcp=2.5
    if not keyword_set(rss) then rss=15.0
    d2r=!DTOR
    nao=!VALUES.F_NAN
    dstep=0.05
  ; dstep=0.02
    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,Rcp,Rss,rr,cth,sth,phrj,Br,Bt,Bp
    brrr=Br & btrr=Bt & bprr=Bp
    sig = -1.
    for is1=0,nstep do begin
        BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
        r0=rr
        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='+-'+strtrim(brrr,2)    ; 26SEP2003
           btrr='+-'+strtrim(btrr,2)    ; 26SEP2003
           bprr='+-'+strtrim(bprr,2)    ; 26SEP2003
           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
        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
           ppb_hc,ga,ha,r,cth,sth,phr,Br,Bt,Bp,nmax=nmax
           bsign=Br/ABS(Br)
           brrr=brrr*bsign
           btrr=btrr*bsign
           bprr=bprr*bsign
           goto, lbl
        endif
        pb0_csss,gc,hc,apar,Rcp,Rss,r,cth,sth,phr,Br,Bt,Bp
    endfor
 lbl:

end