Writes fweight to a HDF5 file
subroutine write_fida_weights
!+ Writes [[libfida:fweight]] to a HDF5 file
!! HDF5 variables
integer(HID_T) :: fid
integer(HSIZE_T), dimension(4) :: dim4
integer(HSIZE_T), dimension(2) :: dim2
integer(HSIZE_T), dimension(1) :: dim1
integer :: error
character(charlim) :: filename
integer :: i,ie,ip,ic,iwav
real(Float64), dimension(:), allocatable :: lambda_arr
real(Float64), dimension(:), allocatable :: ebarr,ptcharr
real(Float64), dimension(:,:), allocatable :: jacobian,e_grid,p_grid
real(Float64), dimension(:,:), allocatable :: vpa_grid,vpe_grid,fida
real(Float64) :: dlambda, wtot, dE, dP
dlambda=(inputs%lambdamax_wght-inputs%lambdamin_wght)/inputs%nlambda_wght
allocate(lambda_arr(inputs%nlambda_wght))
do i=1,inputs%nlambda_wght
lambda_arr(i)=(i-0.5)*dlambda + inputs%lambdamin_wght
enddo
!! define arrays
!! define energy - array
allocate(ebarr(inputs%ne_wght))
do i=1,inputs%ne_wght
ebarr(i)=real(i-0.5)*inputs%emax_wght/real(inputs%ne_wght)
enddo
dE = abs(ebarr(2)-ebarr(1))
!! define pitch - array
allocate(ptcharr(inputs%np_wght))
do i=1,inputs%np_wght
ptcharr(i)=real(i-0.5)*2./real(inputs%np_wght)-1.
enddo
dP = abs(ptcharr(2)-ptcharr(1))
!! define 2d grids
!! define energy grid
allocate(e_grid(inputs%ne_wght,inputs%np_wght))
do i=1,inputs%ne_wght
e_grid(i,:) = ebarr(i)
enddo
!! define pitch grid
allocate(p_grid(inputs%ne_wght,inputs%np_wght))
do i=1,inputs%np_wght
p_grid(:,i) = ptcharr(i)
enddo
!! define velocity space grid
allocate(vpe_grid(inputs%ne_wght,inputs%np_wght)) !! V perpendicular
allocate(vpa_grid(inputs%ne_wght,inputs%np_wght)) !! V parallel
vpa_grid = 100*sqrt((((2.0d3)*e0)/(mass_u*inputs%ab))*e_grid)*p_grid ! [cm/s]
vpe_grid = 100*sqrt((((2.0d3)*e0)/(mass_u*inputs%ab))*e_grid*(1.0-p_grid**2)) ![cm/s]
!! define jacobian to convert between E-p to velocity
allocate(jacobian(inputs%ne_wght,inputs%np_wght))
jacobian = ((inputs%ab*mass_u)/(e0*1.0d3)) *vpe_grid/sqrt(vpa_grid**2 + vpe_grid**2)
!! normalize mean_f
do ic=1,spec_chords%nchan
do ip=1,inputs%np_wght
do ie=1,inputs%ne_wght
wtot = sum(fweight%weight(:,ie,ip,ic))
if((wtot.gt.0.d0)) then
fweight%mean_f(ie,ip,ic) = fweight%mean_f(ie,ip,ic)/wtot
endif
enddo
enddo
enddo
!! Calculate FIDA estimate
allocate(fida(inputs%nlambda_wght,spec_chords%nchan))
do iwav=1,size(fida,1)
fida(iwav,:) = (dE*dP*1d4)*sum(sum(fweight%mean_f(:,:,:)*fweight%weight(iwav,:,:,:),1),1)
enddo
filename=trim(adjustl(inputs%result_dir))//"/"//trim(adjustl(inputs%runid))//"_fida_weights.h5"
!Open HDF5 interface
call h5open_f(error)
!Create file overwriting any existing file
call h5fcreate_f(filename, H5F_ACC_TRUNC_F, fid, error)
dim1(1) = 1
dim2 = [inputs%nlambda_wght, spec_chords%nchan]
dim4 = [inputs%nlambda_wght, inputs%ne_wght, inputs%np_wght, spec_chords%nchan]
call h5ltmake_dataset_int_f(fid,"/nenergy",0,dim1,[inputs%ne_wght], error)
call h5ltmake_dataset_int_f(fid,"/npitch",0,dim1,[inputs%np_wght], error)
call h5ltmake_dataset_int_f(fid,"/nchan",0,dim1,[spec_chords%nchan], error)
call h5ltmake_compressed_dataset_double_f(fid,"/weight",4,dim4,fweight%weight,error)
call h5ltmake_compressed_dataset_double_f(fid,"/fida",2,dim2,fida,error)
call h5ltmake_compressed_dataset_double_f(fid,"/mean_f",3,dim4(2:4),fweight%mean_f,error)
call h5ltmake_compressed_dataset_double_f(fid,"/lambda",1,dim4(1:1),lambda_arr,error)
call h5ltmake_compressed_dataset_double_f(fid,"/energy",1,dim4(2:2),ebarr, error)
call h5ltmake_compressed_dataset_double_f(fid,"/pitch",1,dim4(3:3),ptcharr, error)
call h5ltmake_compressed_dataset_double_f(fid,"/radius",1,dim4(4:4),spec_chords%radius, error)
dim2 = [inputs%ne_wght, inputs%np_wght]
call h5ltmake_compressed_dataset_double_f(fid,"/jacobian",2,dim2, jacobian, error)
call h5ltmake_compressed_dataset_double_f(fid,"/vpe_grid",2,dim2,vpe_grid, error)
call h5ltmake_compressed_dataset_double_f(fid,"/vpa_grid",2,dim2,vpa_grid, error)
call h5ltmake_compressed_dataset_double_f(fid,"/e_grid",2,dim2,e_grid, error)
call h5ltmake_compressed_dataset_double_f(fid,"/p_grid",2,dim2,p_grid, error)
!Add attributes
call h5ltset_attribute_string_f(fid, "/", "version", version, error)
if(inputs%calc_fida_wght.eq.1) then
call h5ltset_attribute_string_f(fid,"/", "description", &
"Line of Sight averaged FIDA E-p space sensitivity/weights " // &
"and spectra calculated by FIDASIM", error)
else
call h5ltset_attribute_string_f(fid,"/", "description", &
"Full FIDA E-p space sensitivity/weights and spectra calculated " // &
"by FIDASIM via Monte Carlo method", error)
endif
call h5ltset_attribute_string_f(fid,"/weight","description", &
"E-p space sensivity/weight of FIDA diagnostic: weight(lambda,energy,pitch,chan)", error)
call h5ltset_attribute_string_f(fid,"/weight","units", &
"(Ph*cm)/(s*nm*sr*fast-ion*dE*dP)",error)
call h5ltset_attribute_string_f(fid,"/fida","units", &
"Ph/(s*nm*sr*m^2)",error )
call h5ltset_attribute_string_f(fid,"/fida","description", &
"Estimate of Fast-ion D-alpha (FIDA) emmision calculated by 1e4*weight*mean_f*dEdP: fida(lambda,chan)", error)
call h5ltset_attribute_string_f(fid,"/mean_f","description", &
"Estimated mean fast-ion distribution function seen by los: mean_f(energy,pitch,chan)", error)
call h5ltset_attribute_string_f(fid,"/mean_f","units", &
"fast-ion/(dE*dP*cm^3)", error)
call h5ltset_attribute_string_f(fid,"/lambda","description", &
"Wavelength array", error)
call h5ltset_attribute_string_f(fid,"/lambda","units","nm", error)
call h5ltset_attribute_string_f(fid,"/nchan", "description", &
"Number of channels", error)
call h5ltset_attribute_string_f(fid,"/nenergy", "description", &
"Number of energy values", error)
call h5ltset_attribute_string_f(fid,"/npitch", "description", &
"Number of pitch value", error)
call h5ltset_attribute_string_f(fid,"/energy","description", &
"Energy array", error)
call h5ltset_attribute_string_f(fid,"/energy", "units","keV", error)
call h5ltset_attribute_string_f(fid,"/pitch", "description", &
"Pitch array: p = v_parallel/v w.r.t. the magnetic field", error)
call h5ltset_attribute_string_f(fid,"/radius", "description", &
"Line of sight radius at midplane or tangency point", error)
call h5ltset_attribute_string_f(fid,"/radius", "units","cm", error)
call h5ltset_attribute_string_f(fid,"/jacobian","description", &
"Jacobian used to convert from E-p space to velocity space", error)
call h5ltset_attribute_string_f(fid,"/jacobian","units", &
"(dE*dP)/(dvpa*dvpe)", error)
call h5ltset_attribute_string_f(fid,"/e_grid","description", &
"2D energy grid", error)
call h5ltset_attribute_string_f(fid,"/e_grid","units","keV", error)
call h5ltset_attribute_string_f(fid,"/p_grid","description", &
"2D pitch grid", error)
call h5ltset_attribute_string_f(fid,"/vpe_grid","description", &
"2D perpendicular velocity grid", error)
call h5ltset_attribute_string_f(fid,"/vpe_grid","units","cm/s", error)
call h5ltset_attribute_string_f(fid,"/vpa_grid","description", &
"2D parallel velocity grid", error)
call h5ltset_attribute_string_f(fid,"/vpa_grid","units","cm/s", error)
!Close file
call h5fclose_f(fid, error)
!Close HDF5 interface
call h5close_f(error)
if(inputs%verbose.ge.1) then
write(*,'(T4,a,a)') 'FIDA weights written to: ', trim(filename)
endif
end subroutine write_fida_weights