1) read with "rfits_n3d.pro"
name='fits/n3d_cr2097P1_m10.fits'
den=rfits_n3d(name, ph,th,r, pbname, index=index)
IDL> help,ph,th,r, pbname, index
PH              DOUBLE    = Array[361]
TH              DOUBLE    = Array[181]
R               DOUBLE    = Array[51]
PBNAME          STRING    = Array[3, 28]
INDEX           STRUCT    = -> <Anonymous> Array[1]


2) read with "readfits.pro" in SSW
n=readfits(name,h)
IDL> help,h,n
H               STRING    = Array[60]
N               LONG      = Array[361, 181, 51]
IDL> print,h
SIMPLE  =                    T / Written by IDL:  Tue Feb 13 19:34:47 2024      
BITPIX  =                   32 /                                                
NAXIS   =                    3 /                                                
NAXIS1  =                  361 /                                                
NAXIS2  =                  181 /                                                
NAXIS3  =                   51 /                                                
CTYPE1  = 'CRLN    '           / Carrington Longitude                           
CRPIX1  =                    1 / first pixel in longitude                       
CRVAL1  =              0.00000 / 0 degree                                       
CDELT1  =              1.00000 / spatial scale in longitude                     
CUNIT1  = 'deg     '           /                                                
CTYPE2  = 'CRLT    '           / Carrington Latitude                            
CRPIX2  =                   91 / Equator                                        
CRVAL2  =              0.00000 / Latitude of Equator                            
CDELT2  =              1.00000 / spatial scale in latitude                      
CUNIT2  = 'deg     '           /                                                
CTYPE3  = 'HECR    '           / Heliocentric distance                          
CRPIX3  =                    1 / first pixel in radial distance                 
CRVAL3  =              1.50000 / 1.5 Rsun                                       
CDELT3  =            0.0500000 / spatial scale in radial distance               
CUNIT3  = 'solRad  '           / solar radius                                   
BTYPE   = 'Density '           / Electron density                               
BUNIT   = 'cm^-3   '           /                                                
INSTRUME= 'SECCHI  '           /                                                
OBSRVTRY= 'STEREO_B'           /                                                
DETECTOR= 'COR1    '           /                                                
DATE-BEG= '2010-05-20T05:55:27.682' / observing time for the first used pb-image
DATE-END= '2010-06-02T12:00:31.213' / observing time for the last used pb-image 
DATE-AVG= '2010-05-26T20:57:59.447' / the middle time for the used pb-images    
DATE    = '2024-02-14T00:34:47.000' / The creation date of the FITS file        
COMMENT The following includes the names of pb-image series used for tomography 
COMMENT /20100520_/055500,055509,055518/_s4c1B.fts                              
COMMENT /20100520_/175500,175509,175518/_s4c1B.fts                              
COMMENT /20100521_/055500,055509,055518/_s4c1B.fts                              
COMMENT /20100521_/175500,175509,175518/_s4c1B.fts                

3. Show the distribution of electron density at a given radial distance (in unit of Rsun)
;to show density map at r=2.5 Rs, density profile along the equator, and radial profile of 
globally averaged density:
n3d_showmap,den,2.5,ph,th,r,nemap,crname='CR2097-P1'
"nemap" return the 2D density at r=2.5 Rs