;+

function bxyz2str, bx, by, bz, $
                   pangle=p, b0angle=b0, $
                   latitude=latin, longitude=cmdin, $
                   radius=radius,pix_size=pix_size, $
                   i_cont=i_contin, flat=flat, point=pointin

;NAME:
;PURPOSE:
;     Put Bx, By, Bz into a field structure
;CATEGORY:
;CALLING SEQUENCE:
;     bstr = bxyz2str(bx,by,bz)
;INPUTS:
;     bx, by, bz: Bx, By, Bz in Gauss
;OPTIONAL INPUT PARAMETERS:
;KEYWORD PARAMETERS
;     i_cont = Continuum image (def = none)
;     pangle = P angle in radians (def = 0)
;     b0angle = B0 angle in radians (def = 0)
;     latitude = latitude in degrees (def = 0)
;     cmd = CMD in degrees (def = 0)
;     radius = solar radius in arcsec (def = 1)
;     pix_size = pixels size in arcsec (def = 1000)
;     point = point structure: take b0angle,pangle,radius,pix_size,cmd,latitude
;             from this structure.
;     /flat = set the lat and cmd arrays to constants.  The default is
;             to set the lat and cmd arrays to the solar values using
;             the pixel size and radius.
;OUTPUTS:
;     bstr = magnetic field structure
;COMMON BLOCKS:
;SIDE EFFECTS:
;RESTRICTIONS:
;PROCEDURE:
;MODIFICATION HISTORY:
;     T. Metcalf 09-Nov-2004
;-


   sz = size(bz)
   if sz[0] NE 2 then message,'ERROR: Bz must have 2 dimensions'
   nx = sz[1]
   ny = sz[2]

   if size(pointin,/tname) EQ 'STRUCT' then begin
      latin    = pointin.lat * !radeg
      cmdin    = pointin.cmd * !radeg
      b0       = pointin.b0
      p        = pointin.p
      pix_size = pointin.pix_size
      radius   = pointin.radius
   endif

   if n_elements(latin) LE 0 then lat = 0.0 else lat = latin
   if n_elements(cmdin) LE 0 then cmd = 0.0 else cmd = cmdin
   if n_elements(p) LE 0 then p = 0.0
   if n_elements(b0) LE 0 then b0 = 0.0
   if n_elements(radius) LE 0 then radius = 1000.
   if n_elements(pix_size) LE 0 then pix_size = [1.,1.]
   if n_elements(pix_size) EQ 1 then pix_size = [pix_size[0],pix_size[0]]

   if n_elements(lat) EQ n_elements(bz) then lat1 = lat[(nx-1)/2,(ny-1)/2] $
   else if n_elements(lat) EQ 1 then lat1 = lat $
   else message,'ERROR: latitude has strange size'

   if n_elements(cmd) EQ n_elements(bz) then cmd1 = cmd[(nx-1)/2,(ny-1)/2] $
   else if n_elements(cmd) EQ 1 then cmd1 = cmd $
   else message,'ERROR: CMD has strange size'

   if n_elements(i_contin) NE n_elements(bz) then i_cont = fltarr(nx,ny) $
   else i_cont = i_contin

   if n_elements(lat) NE n_elements(bz) OR $
      n_elements(cmd) NE n_elements(bz) then begin
      if keyword_set(flat) then begin
         cmd = replicate(cmd1,nx,ny)
         lat = replicate(lat1,nx,ny)
      endif else begin
         ; the arrays are padded to help the interpolation if the 
         ; p angle is not zero.
         xpad = long(nx*(sqrt(2.)-1.0)/2.0+0.5) > 0L
         ypad = long(ny*(sqrt(2.)-1.0)/2.0+0.5) > 0L
         nx1 = nx+2*xpad & ny1 = ny+2*ypad
         xy = lonlat2xy([cmd1,lat1],b0=b0,radius=radius)
         x = pix_size[0]*(lindgen(nx1,ny1) MOD nx1)
         y = pix_size[1]*(lindgen(nx1,ny1) / nx1)
         x = x - x[(nx1-1)/2,(ny1-1)/2] + xy[0]
         y = y - y[(nx1-1)/2,(ny1-1)/2] + xy[1]
         ; rotate x and y by the p angle then compute the lat and lon arrays
         if keyword_set(p) then begin
            x = rot(x,-p*!radeg,1.0,(nx1-1)/2,(ny1-1)/2,cubic=-0.5,/pivot)
            y = rot(y,-p*!radeg,1.0,(nx1-1)/2,(ny1-1)/2,cubic=-0.5,/pivot)
         endif
         ll = xy2lonlat([[x],[y]],b0=b0,radius=radius)
         ll = xy2lonlat(transpose(reform([reform(x,nx1*ny1),reform(y,nx1*ny1)],nx1*ny1,2)), $
                        b0=b0,radius=radius)
         cmd = (reform(ll[0,*],nx1,ny1))[xpad:xpad+nx-1,ypad:ypad+ny-1]
         lat = (reform(ll[1,*],nx1,ny1))[xpad:xpad+nx-1,ypad:ypad+ny-1]
      endelse
   endif

   bxyz,bxi,byi,bzi,bx,by,bz,p=p,b0=b0,lat=lat,cmd=cmd,/helio2image
   bl = bzi
   bt = sqrt(bxi^2+byi^2)
   ba = atan(-bxi,byi)*!radeg

   stime = !stime
   pname = strmid(stime,0,2)+strmid(stime,3,3)+strmid(stime,9,2)+ $
           strmid(stime,12,2)+strmid(stime,15,2)+strmid(stime,18,2)
   if strmid(pname,0,1) EQ ' ' then pname = '0'+strmid(pname,1)
   pname = strcompress('P'+pname,/remove_all)  ; structure name
   point = create_struct(name=pname,['b0','p','cmd','lat','pix_size','radius'], $
                                    b0,p,cmd1*!dtor,lat1*!dtor,pix_size,radius)

   b = {b_long: bl, b_trans: bt, b_azim: ba, $
        bx:bx, by:by, bz:bz, $
        latitude:lat, cmd: cmd, $
        i_cont:i_cont,point: point}

   return,jxyz(b,/quiet)


end