;
; 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,'SYNTAX - pro gh_mbrsc,mbrsc,Nmax,gm,hm,ghmf=,/ifl'
       print,'INPUT  - mbfsc,Nmax,ghmf,/ifl'
       print,'OUTPUT - ghmf,gm,hm'
       print,'NOTE   - mbfsc is an array or, if /ifl, a filename' 
       return
    endif

    if keyword_set(ifl) then begin
       bra=readfits(mbrsc,header) 
       tt=mbrsc
    ;  bra=congrid(bra,360,180)
    endif else begin
       tt=' '
       bra=mbrsc
    endelse
    sz=SIZE(bra)
    xsz=sz(1) & ysz=sz(2)
    zgrid,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,tt
       printf,unit,Nmax
    endif
    gm=fltarr(Nmax+1,Nmax+1)
    hm=gm
    ij=xsz*ysz
    for n=1,Nmax do begin         ; n=1 because no monopole at photosphere
        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