;
; PURPOSE:   Compute GENERAL G & H based on WSO data
;              wbrsc (Br array or WSC filename) and Nmax 
; WRITTEN:   03May2000  Xuepu
;
pro pgh_wbrsc,wbrsc,Nmax,gw,hw,GHWF=ghwf,IFL=ifl,CALI=cali

    if N_params() LT 1 then begin
       print,'SYNTAX - pro gh_wbrsc,wbrsc,Nmax,ghwf,gw,hw,/ifl,/cali'
       print,'INPUT  - wbrsc,Nmax,ghwf,/ifl,/cali'
       print,'OUTPUT - ghwf,gw,hw'
       print,'NOTE   - wbrsc is an array or, if /ifl, a filename' 
       return
    endif

    if keyword_set(ifl) then begin
       tt=wbrsc
       if keyword_set(cali) then bra=readwsosc(wbrsc,/cali,/br)*0.01 $
                            else bra=readwsosc(wbrsc,/br)*0.01 
    endif else begin
       tt=' '
       bra=wbrsc
    endelse
    zgrid,72,30,phd,thd,lad,cth,sth
    pdp,cth,Nmax,P,dP
    csmph,phd,nmax,cmph,smph

    if keyword_set(ghwf) then begin
       get_lun,unit
       openw,unit,ghwf
       printf,unit,tt
       printf,unit,Nmax
    endif
    gw=fltarr(Nmax+1,Nmax+1)
    hw=gw
    ij=72*30
    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
            gw(n,m)=0
            hw(n,m)=0
            for i=0,71 do begin
                for j=0,29 do begin
                    gw(n,m)=gw(n,m)+bra(i,j)*P(n,m,j)*cmph(m,i)
                    hw(n,m)=hw(n,m)+bra(i,j)*P(n,m,j)*smph(m,i)
                endfor
            endfor
            gw(n,m)=coe*gw(n,m)
            hw(n,m)=coe*hw(n,m)
            if keyword_set(ghwf) then printf,unit,n,m,gw(n,m),hw(n,m)
        endfor
    endfor
    if keyword_set(ghwf) then free_lun,unit

end