;
; 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,phri,thrj,brrr,btrr,bprr,APAR=apar,RCP=rcp,RSS=rss

; print,'This is /home/xuepu/ZIDLLIB/CSSS/ppb_csss.pro'
    if N_params() lt 1 then begin
       print,'ppb_csss,ga,ha,gc,hc,rr,phrj,thri,brrr,btrr,bprr,apar=,rcp=,rss='
       print,'OUTPUT: brrr,btrr,bprr'
       return
    endif
    if not keyword_set(apar) then apar=0.2
    if not keyword_set(rcp) then rcp=2.25
    if not keyword_set(rss) then rss=15.0
    d2r=!DTOR
    nao=!VALUES.F_NAN
    dstep=0.01
    nstep=500
; ***
    r=rr
    cth=cos(thrj)
    sth=sin(thrj)
    cph=cos(phri)
    sph=sin(phri)
    x = r*sth*cph
    y = r*sth*sph
    z = r*cth
    pb0_csss,gc,hc,rr,phri,cth,sth,apar,Rcp,Rss,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
           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,phr,cth,sth,Br,Bt,Bp
           bsign=Br/ABS(Br)
           brrr=brrr*bsign
           btrr=btrr*bsign
           bprr=bprr*bsign
           goto, lbl
        endif
        pb0_csss,gc,hc,r,phr,cth,sth,apar,Rcp,Rss,Br,Bt,Bp
    endfor
 lbl:

end