;+  Convert fits header for n3d-fits file into header structure
;               
; Use         : IDL> index = FITSHD2INDEX(fits_hdr)
;    
; Inputs      : fits_hdr	FITS header, STRARR
;               
; Outputs     : index   	header structure.
;               
; Calls       : FXPAR, GETTOK
;
; Category    : Data_Handling, I_O
;
; History:
;  Adapted  to add a keyword "name" to read names of pB-series in the fits-header             
;-            

;__________________________________________________________________________________________________________
;

FUNCTION FITSHD2INDEX, fts_hdr, name=name

   len = N_ELEMENTS(fts_hdr)

      key = STRTRIM(GETTOK(GETTOK(fts_hdr(0), ' '),'='),2)
      val = FXPAR(fts_hdr, key)
      hdr = CREATE_STRUCT(key, val)
      FOR i=1, len-1 DO BEGIN  	;** for each keyword in the fits header
         key = STRTRIM(GETTOK(GETTOK(fts_hdr(i), ' '),'='),2)help,f
                                	;** get the keyword
         val = FXPAR(fts_hdr, key)	;** and get its value
         IF (DATATYPE(val(0)) EQ 'STR') THEN val = STRTRIM(val(0), 2)
         IF (key EQ 'DATE-BEG') THEN key = 'DATE_BEG'
         IF (key EQ 'DATE-END') THEN key = 'DATE_END'
	 IF (key EQ 'DATE-AVG') THEN key = 'DATE_AVG'
         IF (key NE 'END') and (key NE '' and key NE 'COMMENT') THEN $
            hdr = CREATE_STRUCT(hdr, key, val)	
              ;** and add to the structure definition
      ENDFOR
       
      if keyword_set(name) then begin
	 tmp_hdr=fts_hdr
 	 w =where(GETTOK(tmp_hdr, ' ') eq 'COMMENT', nw)
         cmt =fts_hdr[w]
	 name=strarr(3,nw-1)
	 for i=1,nw-1 do begin
   	     a=str_sep(cmt[i],'/')
	     b=str_sep(a[2],',')
	     name[0,i-1]=a[1]+b[0]+a[3]
	     name[1,i-1]=a[1]+b[1]+a[3]
	     name[2,i-1]=a[1]+b[2]+a[3]
         endfor
	 name=strcompress(name,/remove_all)
	 hdr = CREATE_STRUCT(hdr, 'filename', name)
      endif
    ; hdr_tags = TAG_NAMES(hdr)

   RETURN, hdr

END


 FUNCTION rfits_n3d, filename, ph,th,r, pbname, index=index
 ;+ read n3d fits file and return 3D cubic density data in cm^-3
 ; Input: 
 ;  filename: standard fits file
 ; Return:
 ;   data=rfits_n3d(filename, ph,th,r)
 ;  data[m,n,k] -3D cube with dimensions for [lon, lat, r]
 ; Output:
 ;  ph[m] Carrington longitude in deg
 ;  th[n] Carrington latitude in deg
 ;  r[n] radial distance in Rsun
 ;   pbname[3,*], output in string array, including names for pb image sereis
 ;     used for the density reconstruction by tomography
 ;  Keyword:
 ;   index - fits header in index struct
 ;-
  data=readfits(filename,hdr)
  if n_params() gt 4 then begin
     print,'return filenames for pB images'
     index=fitshd2index(hdr,/name) 
     pbname=index.filename
  endif $   
  else index=fitshd2index(hdr)

  ss=size(data)
  iph=findgen(ss[1])+1
  ith=findgen(ss[2])+1
  ir=findgen(ss[3])+1
  
  ph=(iph-index.crpix1)*index.cdelt1+index.crval1
  th=(ith-index.crpix2)*index.cdelt2+index.crval2
  r=(ir-index.crpix3)*index.cdelt3+index.crval3

  return,double(data)
 END