;
; PURPOSE:   Compute GENERAL G & H based on MDI or KPK Br SC data,
;              mbrsc (Br array or MDI filename) and Nmax 
; WRITTEN:   09May2000  Xuepu
;
pro pgh_mbrsc,mbrsc,nmax,gm,hm,GHMF=ghmf,IFL=ifl

    if N_params() LT 1 then begin
       print,'pro pgh_mbrsc,mbrsc,nmax,gm,hm,ghmf=,/ifl'
       print,'INPUT  - mbfsc,nmax,ghmf=,/ifl'
       print,'OUTPUT - gm,hm,ghmf'
       print,'NOTE   - mbfsc is an array or a filename when /ifl'
       return
    endif

    if keyword_set(ifl) then begin
       bra=readfits(mbrsc,header) 
    endif else begin
       bra=mbrsc
    endelse
    sz=SIZE(bra)
    xsz=sz(1) & ysz=sz(2)
    pzgrid,xsz,ysz,phd,thd,lad,cth,sth
    pdp,cth,nmax,P,dP
    csmph,phd,nmax,cmph,smph
    if keyword_set(ghmf) then begin
       get_lun,unit
       openw,unit,ghmf
       printf,unit,nmax
    endif
    gm=fltarr(nmax+1,nmax+1)
    hm=gm
    ij=xsz*ysz
    for n=0,nmax do begin         ; n=0 corresponding to a monopole 
        coe=(2*n+1.0)/ij
        for m=0,n do begin
            gm(n,m)=0
            hm(n,m)=0
            for i=0,xsz-1 do begin
                for j=0,ysz-1 do begin
                    gm(n,m)=gm(n,m)+bra(i,j)*P(n,m,j)*cmph(m,i)
                    hm(n,m)=hm(n,m)+bra(i,j)*P(n,m,j)*smph(m,i)
                endfor
            endfor
            gm(n,m)=coe*gm(n,m)
            hm(n,m)=coe*hm(n,m)
            if keyword_set(ghmf) then printf,unit,n,m,gm(n,m),hm(n,m)
        endfor
    endfor
    if keyword_set(ghmf) then free_lun,unit

end