;
; MODIFIED:	Xuepu Zhao	13MAR2006
;
pro psb_hc,ga,ha,r,sth,cmph,smph,P,dP,bbr,bbt,bbp,APAR=apar

    if N_params() LT 1 then begin
       print,'pro sb_hc,r,sth,cmph,smph,P,dP,bbr,bbt,bbp,apar='
       return
    endif
    if not keyword_set(apar) then apar=0.2
    sz=SIZE(cmph) & nmax=sz(1)-1 & iph=sz(2)
    jth=N_ELEMENTS(sth)
    Gfi = fltarr(nmax+1,nmax+1)
    dGfi = Gfi
    pj = Gfi
    dpj = pj
    bbr = fltarr(iph,jth)
    bbt = bbr
    bbp = bbr
    hcHf,nmax,apar,r,Hr,dHr,d2Hr
    Gfn,ga,ha,cmph,smph,Gf,dGf
    FOR i = 0, iph-1 do begin
        Gfi(*,*) = Gf(*,*,i)
        dGfi(*,*) = dGf(*,*,i)
        FOR j = 0, jth-1 do begin
            sthj = sth(j)
            pj(*,*) = P(*,*,j)
            dpj(*,*) = dP(*,*,j)
            brtp,nmax,apar,r,sthj,Hr,dHr,pj,dpj,Gfi,dGfi,Br,Bt,Bp
            bbr(i,j) = Br
            bbt(i,j) = Bt
            bbp(i,j) = Bp
        ENDFOR
    ENDFOR

end