;
; PURPOSE: Calculate BrBtBp on a spherical surface located at r
; INPUT:   g,h,gc,hc,r,phd,thd,nmax,apar,Rcp,Rss
; OUTPUT:  bbr,bbt,bbp
; MODIFIED: 08DEC2005	Xuepu Zhao
;
pro psb_csss,g,h,gc,hc,rr,phd,thd,nmax,apar,Rcp,Rss,bbr,bbt,bbp
  
    if N_params() lt 1 then begin
       print,'psb_csss,g,h,gc,hc,rr,phd,thd,nmax,apar,Rcp,Rss,bbr,bbt,bbp'
       return
    endif

    ix = N_ELEMENTS(phd)
    jy = N_ELEMENTS(thd)
    bbr = fltarr(ix,jy)
    bbt = bbr
    bbp = bbr
    phr = phd*!DTOR
    thr = thd*!DTOR
    FOR i = 0, ix-1 do begin
        phri=phr(i)
        FOR j = 0, jy-1 do begin
            thrj=thr(j)
            ppb_csss,g,h,gc,hc,rr,thrj,phri,nmax=9,br,bt,bp
            bbr(i,j) = br
            bbt(i,j) = bt
            bbp(i,j) = bp
        ENDFOR
    ENDFOR

end