Loading sources/fit_data00.f90deleted 100644 → 0 +0 −584 Original line number Diff line number Diff line !! ========================================================================== !! module fit_data !! ========================================================================== module fit_data use constants_module use data_types use echo_shapes use dump_data use fit_utils use mathutil use logger implicit none !!mm>> new to save offset info for reporting (experimental) integer, parameter :: MAX_TOTAL_POINTS = 10000 real(kind=DBL) :: DJOFFSET_LIMIT = 0.2d-6 integer :: count_dj_offset = 0 integer :: runno_point_stat(3,MAX_TOTAL_POINTS) = 0 real(kind=DBL) :: dj_phaseoffsets(MAX_TOTAL_POINTS) = 0 !!mm<< contains !! ======================================================================== !> @brief Fit echo data !! !! @param [in,out] data - scan_struct data (a @f$ \tau @f$ scan) !! @param [in] fit_flag - fitting level flag @see "phase_" flags in drspine_parameters !! @param [in,optional] phase_offset - initial phase offset !! !! @description !! This is user callable subroutine - public interface to fit_data module. subroutine fit_echo_data(data, fit_flag, phase_offset) type(scan_struct), target, intent(inout) :: data integer, intent(in) :: fit_flag real(kind=DBL), intent(in), optional :: phase_offset ! type(scan_data_struct), pointer :: sig_point, ref_point integer :: iscan real(kind=DBL) :: djoffset call msg_info('fit_echo_data', 'fitting echo data '//trim(data%file)//& ', role='//trim(cformat_role(data%role))//& ', flag='//trim(cformat_fit_flag(fit_flag))) djoffset = 0.0d0 if ( present(phase_offset) .and. data%number_of_points>0 ) & djoffset = get_field_integral(deg2rad(phase_offset), data%scan_point(1)%physics%lambda) !just need a wavelength ref_point => null() do iscan=1, data%number_of_points call msg_debug('fit_echo_data', msg_fmt("('scan point[0]: ', i4)", iscan)) sig_point => data%scan_point(iscan) if ( associated(sig_point%matching_ref) ) ref_point => sig_point%matching_ref !! call fit_center_bin(sig_point, ref_point, fit_flag) if ( iand(fit_flag, PHASE_FIT) /= 0 ) then ! fit phase for full b-w (no smoothing yet) call fit_echo_param(sig_point, ref_point, iand(fit_flag, not(PHASE_FIT_SMOOTH))) ! smooth if ( iand(fit_flag, PHASE_FIT_SMOOTH) /= 0 ) then call fit_echo_param(sig_point, sig_point, ior(PHASE_FIT, PHASE_FIT_SMOOTH)) endif if ( .not. associated(sig_point%matching_ref) ) sig_point%matching_ref => sig_point ref_point => sig_point endif if ( associated(sig_point%matching_ref) ) then call get_echo_param(sig_point, sig_point%matching_ref, .false.) ! extract if ( iand(fit_flag, PHASE_FIT_OFFSET) /=0 ) then if (iand(fit_flag, PHASE_FIT_RESET) /=0 ) djoffset=0 call get_phase_offset(sig_point, sig_point%matching_ref, djoffset) call get_echo_param(sig_point, sig_point%matching_ref, .true., djoffset) end if end if end do end subroutine fit_echo_data !! ======================================================================== !> @brief This subroutine fits phase scans (for resolution samples) !! !! @param dat_point - scan_data_struct pointer with phase scan data !! @param ref_point - scan_data_struct pointer with reference data flag (may be null) !! @param fit_flag - fitting level flag @see "phase_" flags in drspine_parameters !! !! @description !! This subroutine will perform symmetry phase scan fit over entire scan_data_struct !! (typically all detector pixels). !! It loops from the center of the detector in a spiral sequence. !! !! The supported flags are !! - PHASE_FIT - (implied) will perform phase fitting !! - PHASE_FIT_HOPPING - fit phase (full wavelength band) with "hopping" !! - PHASE_FIT_AMPPOS - fitted amplitude must be positive !! - PHASE_FIT_AMPNEG - fitted amplitude must be negative !! - PHASE_FIT_SMOOTH - smooth phase map and refit !! !! @note If ref_point is not null (i.e. associated) it will be used as a starting point !! for fitting (unless PHASE_FIT_SMOOTH is set). !! @note The ref_point may refer to the data itself. subroutine fit_echo_param(dat_point, ref_point, fit_flag) type(scan_data_struct),pointer :: dat_point, ref_point integer, intent(in) :: fit_flag ! !! it=0 fit phases only for the entire wavelength band integer, parameter :: it=0 ! fit full b-w integer :: ix, iy, ipix integer :: xpix, ypix type(lambda_struct) :: eshape_lams type(phase_scan_struct),pointer :: datpix, refpix real(kind=DBL) :: dj0 ! symmetry phase initial guess real(kind=DBL) :: ddj ! symmetry phase limits dj0 +/- ddj (typically +/- 90 degrees in phase) real(kind=DBL) :: cdj ! smoothing correction unit typically corresponding to 360 degree in phase ! integer, dimension(:), allocatable :: ixav integer, dimension(:), allocatable :: iyav !call print_image(97, scan_point%pixelbin(it,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAP1', z6)", fit_flag), .FALSE.) ddj = get_field_integral(reduction_parameters%symmetry_phase_max_dev, dat_point%physics%lambda) cdj = get_field_integral(reduction_parameters%symmetry_phase_corr_step, dat_point%physics%lambda) if ( iand(fit_flag, PHASE_FIT_SMOOTH)/=0 ) then call smooth_phase_map(dat_point%pixelbin(it,:,:)%delta_J_symm, & reduction_parameters%symmetry_phase_patch_size, ddj, cdj) !call print_image(97, dat_point%pixelbin(it,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAPS', z7)", fit_flag), .FALSE.) !return endif xpix = max(1, dat_point%no_xpix) ! at least one ypix = max(1, dat_point%no_ypix) allocate(ixav(xpix*ypix)) allocate(iyav(xpix*ypix)) call create_spiralad(ixav,iyav, xpix, ypix) call init_lambda_struct(eshape_lams, dat_point%spectrum(0)) pixel_loop: do ipix=1, xpix*ypix dj0 = 0.0_DBL refpix => null() !if (ipix==0) then ! datpix => dat_point%centerbin ! if ( associated(ref_point) ) refpix => ref_point%centerbin !else ix = ixav(ipix) iy = iyav(ipix) datpix => dat_point%pixelbin(it, ix, iy) if ( associated(ref_point) ) refpix => ref_point%pixelbin(it, ix, iy) !end if if ( iand(fit_flag, PHASE_FIT_SMOOTH)/=0 .and. ipix/=0 ) then dj0 = datpix%delta_J_symm%value else if ( associated(refpix) ) then dj0 = refpix%delta_J_symm%value endif call fit_phase_scan(datpix, eshape_lams, fit_flag, dj0, ddj) call msg_debug('fit_phase', clog_fit_params(datpix, 'pha:')) end do pixel_loop deallocate(ixav) deallocate(iyav) !call print_image(97, scan_point%pixelbin(it,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAP2', z6)", fit_flag), .FALSE.) end subroutine fit_echo_param !! ======================================================================== !> @brief This subroutine fits phase scans (for resolution samples) !! !! @param dat_point - scan_data_struct pointer with phase scan data !! @param ref_point - scan_data_struct pointer with reference data flag (may be null) !! @param fit_flag - fitting level flag @see "phase_" flags in drspine_parameters !! subroutine fit_center_bin(dat_point, ref_point, fit_flag) type(scan_data_struct),pointer :: dat_point, ref_point integer, intent(in) :: fit_flag ! type(lambda_struct) :: eshape_lams type(phase_scan_struct),pointer :: datpix, refpix real(kind=DBL) :: dj0 ! symmetry phase initial guess real(kind=DBL) :: ddj ! symmetry phase limits dj0 +/- ddj (typically +/- 90 degrees in phase) ! ddj = get_field_integral(reduction_parameters%symmetry_phase_max_dev, dat_point%physics%lambda) call init_lambda_struct(eshape_lams, dat_point%spectrum(0)) dj0 = 0.0_DBL refpix => null() datpix => dat_point%centerbin if ( associated(ref_point) ) then refpix => ref_point%centerbin dj0 = refpix%delta_J_symm%value end if call fit_phase_scan(datpix, eshape_lams, fit_flag, dj0, ddj) call msg_debug('fit_phase', clog_fit_params(datpix, 'pha:')) end subroutine fit_center_bin !! ============================================================================ !> @brief This subroutine gets (extracts) echo parameters using ref_point phase !! !! @param scan_point - scan_data_struct pointer with phase scan data !! @param ref_point - scan_data_struct pointer with reference data flag (cannot be null) !! !! @description !! !! @note The ref_point may refer to the data itself. subroutine get_echo_param(scan_point, res_point, apply_offset, offset) type(scan_data_struct) :: scan_point, res_point logical, intent(in) :: apply_offset real(kind=DBL), intent(in), optional :: offset ! integer :: it, ix, iy integer :: tpix, xpix, ypix type(lambda_struct) :: eshape_lams type(value_struct) :: dj real(kind=DBL) :: djoffset tpix = scan_point%no_lam xpix = scan_point%no_xpix ypix = scan_point%no_ypix call init_lambda_struct(eshape_lams, scan_point%spectrum(0)) djoffset = 0 if (apply_offset) then if (present(offset) ) djoffset = offset call msg_info('get_echo_param', & trim(msg_fmt("(1x,'applying phase offset: ',g12.5)", djoffset))//& trim(msg_fmt("(1x,'to scan[',i0,'/',i0,']:')", [scan_point%id, scan_point%point]))) !!mm>> new to save offset info for reporting (experimental) if(count_dj_offset < MAX_TOTAL_POINTS) count_dj_offset = count_dj_offset + 1 runno_point_stat(1,count_dj_offset) = scan_point%id runno_point_stat(2,count_dj_offset) = scan_point%point runno_point_stat(3,count_dj_offset) = int(abs(djoffset)/DJOFFSET_LIMIT) dj_phaseoffsets ( count_dj_offset) = djoffset !!mm<< end if do it=0, tpix !call msg_debug('get_echo_param', & ! msg_fmt("('scan point[',i2,']: ', i4)", [it ,iscan])) do ix=1, xpix do iy=1, ypix ! FIXME what if res and scan are one and the same? if ( res_point%pixelbin( 0, ix, iy)%status .ne. PIXEL_OK ) then scan_point%pixelbin(it, ix, iy)%status = & ior(scan_point%pixelbin(it, ix, iy)%status, & PIXEL_NO_RESOLUTION) cycle ! FIXME: shall we allow to fit without resolution?? endif dj = res_point%pixelbin(0, ix, iy)%delta_J_symm !if (apply_offset) then ! dj = add_value(dj, off) !end if !dj0 = res_point%pixelbin(0, ix, iy)%delta_J_symm%value !dje = sqrt(res_point%pixelbin(0, ix, iy)%delta_J_symm%sigma2) !call fit_phase_scan(scan_point%pixelbin(it, ix, iy), eshape_lams, PHASE_USE, dj0, dje) call fit_phase_scan(scan_point%pixelbin(it, ix, iy), eshape_lams, PHASE_USE, dj%value+djoffset, sqrt(dj%sigma2)) call msg_debug('get_echo_param', clog_fit_params(scan_point%pixelbin(it, ix, iy), 'ext:')) end do end do end do !call print_image(97, scan_point%pixelbin(0,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAP', i4)", iscan), .FALSE.) end subroutine get_echo_param !! ======================================================================== !> @brief Perform actual symmetry phase fit !! !! @param phase_scan - phase scan data !! @param eshape_lams - wavelength spectrum !! @param fit_flag - fit level flag !! @param dj0 - initial guess for the symmetry phase !! @param djlim - phase minimization limits in Tm (Tesla-meters). !! @param tolerance - (optional) minimization tolerance !! !! @description !! Algorithm outline: !! - pre-fit check (only if fitting and not smoothing) !! - check if we have enough "wiggling" in the echo !! - check if the average in inside up/down bounds !! - fit the phase withing the limits (if PHASE_FIT is set) !! - extract echo parameters (average, amplitude, phase) with errors !! - post-fit check !! - check @f$\chi^2@f$ !! - check fit amplitude and average !! subroutine fit_phase_scan(phase_scan, eshape_lams, fit_flag, dj0, dj0_lim, tolerance) type(phase_scan_struct), intent(inout) :: phase_scan ! input data type(lambda_struct), intent(inout) :: eshape_lams ! wavelength spectrum integer, intent(in) :: fit_flag ! fit flag real(kind=DBL), intent(in) :: dj0 ! initial guess real(kind=DBL), intent(in) :: dj0_lim ! fit limits real(kind=DBL), intent(in), optional :: tolerance ! fit tolerance ! integer, parameter :: NHOPS = 3 real(kind=DBL), dimension(-NHOPS:NHOPS,2) :: dj ! dj bracket(s) real(kind=DBL) :: dj0_fit, dj0_err ! fitted phase and error real(kind=DBL) :: dj0_tol ! tolerance real(kind=DBL) :: chi2_fit ! absolute chi2 minimum real(kind=DBL) :: xchi, xdj ! tmp vars integer :: stat ! minimization status integer :: i real(kind=DBL), dimension(2) :: dj0_hard ! type(value_struct) :: amp ! optional parameters dj0_tol = reduction_parameters%symmetry_phase_tolerance if ( present(tolerance) ) dj0_tol = tolerance if(phase_scan%amplitude%state == VAL_UNDEFINED) then ! (PAZ: FIXME, weird way to test for "not checked" raw data) if ( iand(fit_flag, PHASE_FIT) > 0 .and. iand(fit_flag, PHASE_FIT_SMOOTH)==0 ) then call check_raw_data(phase_scan, reduction_parameters%min_sigma_relative) end if endif if (phase_scan%status/=PIXEL_OK) return eshape_lams%nbin = [phase_scan%tbin_1, phase_scan%tbin_2] dj0_fit = dj0 dj0_err = HUGE(dj0) chi2_fit = HUGE(dj0) dj0_hard(1) = MINVAL(phase_scan%delta_J(1:phase_scan%n_down)%value) dj0_hard(2) = MAXVAL(phase_scan%delta_J(1:phase_scan%n_down)%value) if ( iand(fit_flag, PHASE_FIT) /= 0 ) then if ( iand(fit_flag, PHASE_FIT_HOPPING) == 0 ) then dj(0,1) = MAX(dj0_hard(1), dj0 - dj0_lim) ! brent needs a bracket dj(0,2) = MIN(dj0_hard(2), dj0 + dj0_lim) if ( dj(0,1)>dj(0,2) ) call swap(dj(0,1), dj(0,2)) stat = linextpha1(eshape, phase_scan, eshape_lams, dj(0,:), dj0_tol) !minimize if (stat<0) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED) return end if dj0_fit = 0.5*(dj(0,1)+dj(0,2)) else do i=-NHOPS,NHOPS dj(i,1) = MAX(dj0_hard(1), 2*(i-1)*dj0_lim) ! brent needs a bracket dj(i,2) = MIN(dj0_hard(2), 2*(i+1)*dj0_lim) end do ! fit with 'hopping' stat = linextpha2(eshape, phase_scan, eshape_lams, dj, dj0_tol) !minimize if (stat<0) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED) return end if hops: do i=-NHOPS,NHOPS xdj = 0.5*(dj(i,1)+dj(i,2)) call linextract(eshape, phase_scan, eshape_lams, xdj, xchi, amp=amp) if (iand(fit_flag, PHASE_FIT_AMPPOS)/=0 .and. amp%value<0) cycle hops if (iand(fit_flag, PHASE_FIT_AMPNEG)/=0 .and. amp%value>0) cycle hops if ( xchi<chi2_fit ) then dj0_fit = xdj chi2_fit = xchi end if end do hops end if end if if (phase_scan%status/=PIXEL_OK) return call get_phase_fit_values(eshape, phase_scan, eshape_lams, dj0_fit) ! get fit values phase_scan%fqt = get_fqt(phase_scan) call check_fit_par(phase_scan, fit_flag, & reduction_parameters%fqt_min_value,& reduction_parameters%average_max_deviation, & reduction_parameters%max_chisquare) end subroutine fit_phase_scan ! ======================================================================== !> Helper subroutine. !! Check if we have enough "wiggling" in the echo signal by checking !! simga/average of the data against a threshold !! Use this function *before* fitting phase !! @param phase_scan - phase_scan_struct data !! @param threshold - min value for sigma/average subroutine check_raw_data(phase_scan, threshold) type(phase_scan_struct), intent(inout) :: phase_scan real(kind=DBL), intent(in) :: threshold ! real(kind=DBL) :: ave, sig, up, dn ! ave = phase_scan%average%value sig = sqrt(phase_scan%average%sigma2) !! up = max(phase_scan%up%value, phase_scan%down%value) !! dn = min(phase_scan%up%value, phase_scan%down%value) up = phase_scan%up%value dn = phase_scan%up%value ! check counting statistics if ( phase_scan%average_raw%value.lt.reduction_parameters%min_counts_per_pixel) then phase_scan%status = ior(phase_scan%status, PIXEL_STATISTICS) !call msg_debug('check_raw_data', trim(phase_scan%cpixel)//' not enough statistics'//& ! trim(msg_fmt('(" counts=",g12.6)', phase_scan%average_raw%value))//& ! trim(msg_fmt('(" status=",i0 )', phase_scan%status ))) return end if ! if (ave<0) then phase_scan%status = ior(phase_scan%status, PIXEL_AMPLITUDE_RAW) else if (sig/ave<threshold) & phase_scan%status = ior(phase_scan%status, PIXEL_AMPLITUDE_RAW) if (ave<dn .or. up<ave) & phase_scan%status = ior(phase_scan%status, PIXEL_AVERAGE_OUT_OF_BOUNDS) end if end subroutine check_raw_data ! ======================================================================== !> Helper subroutine. !! 1. Check if the relative amplitude of the fit (udratio) is high enough !! 2. Check if the fit average is consistent with up/down average !! !! Use this function *after* fitting phase !! @param phase_scan - phase_scan_struct data !! @param fit_flag - fitting flag !! @param threshold_amp - min value for sigma/average !! @param threshold_ave - max relative deviation !! @param threshold_chi2 - max allowed chi2 subroutine check_fit_par(phase_scan, fit_flag, threshold_amp, threshold_ave, threshold_chi2) type(phase_scan_struct), intent(inout) :: phase_scan integer, intent(in) :: fit_flag real(kind=DBL), intent(in) :: threshold_amp real(kind=DBL), intent(in) :: threshold_ave real(kind=DBL), intent(in) :: threshold_chi2 ! real(kind=DBL) :: ave, amp, fqt, uda ave = phase_scan%average%value amp = phase_scan%amplitude%value fqt = phase_scan%fqt%value uda = (phase_scan%up%value + phase_scan%down%value)/2 if ( phase_scan%chi_squared > threshold_chi2 ) & phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_CHISQUARED) if ( iand(fit_flag, PHASE_FIT_AMPPOS) > 0 .and. amp<0 ) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_NEGATIVE) endif if ( iand(fit_flag, PHASE_FIT_AMPNEG) > 0 .and. amp>0 ) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_POSITIVE) endif if ( iand(fit_flag, PHASE_FIT) > 0 ) then if ( ave>0 .and. abs(fqt)<threshold_amp) & phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_AMPLITUDE) if ( abs(ave-uda)/uda > threshold_ave ) & phase_scan%status = ior(phase_scan%status, PIXEL_AVERAGE_INCONSISTENT) endif return end subroutine check_fit_par ! ======================================================================== !> Helper subroutine. !! Log fit parameters function clog_fit_params(phase_scan, cprefix) result(cline) character(len=MAX_LINE_LENGTH) :: cline type(phase_scan_struct), intent(in) :: phase_scan character(len=*), intent(in) :: cprefix write(cline,'(a,1x,a,2(1x,g13.6),4(1x,a),1x,g13.6,z8)') & trim(cprefix), & trim(phase_scan%cpixel), & !indices(1), indices(2), & phase_scan%Q_abs*ANGSTROEM , & phase_scan%tau/NS, & trim(cformat_value(phase_scan%average, 'g13.6', MAX_NAME_LENGTH)),& trim(cformat_value(phase_scan%amplitude,'g13.6', MAX_NAME_LENGTH)),& trim(cformat_value(phase_scan%delta_J_symm, 'g13.6', MAX_NAME_LENGTH)),& trim(cformat_value(phase_scan%fqt, 'g13.6', MAX_NAME_LENGTH)),& phase_scan%chi_squared, & phase_scan%status return end function clog_fit_params function select_best_phase(sfun, phase_scan, eshape_lams, dJ0) result(res) procedure(echo_shape_function) :: sfun type(phase_scan_struct), intent(in) :: phase_scan type(lambda_struct), intent(in) :: eshape_lams real(kind=DBL), intent(in), dimension(:) :: dJ0 real(kind=DBL) :: res ! real(kind=DBL) :: chi2, chi2_min type(value_struct) :: amp integer :: k ! res = 0.0_DBL chi2_min = HUGE(1.0_DBL) !print *, phase_scan%delta_J(1)%value, phase_scan%delta_J(phase_scan%n_down)%value do k=1, SIZE(dJ0) call linextract(sfun, phase_scan, eshape_lams, dJ0(k), chi2, amp=amp) !print *, k, dJ0(k), chi2, amp%value if ( dJ0(k) < phase_scan%delta_J(1)%value .or. phase_scan%delta_J(phase_scan%n_down)%value < dJ0(k)) cycle if ( amp%value<0) cycle !print *, 'OK' if ( chi2 < chi2_min ) then res = dJ0(k) chi2_min = chi2 end if end do end function select_best_phase !!mm>> !! ======================================================================== !!> @brief Check_Dj_Offset(run,ntau, dj_variance, dj_value) !! !! @param [in] run - numor of run to be considered !! @param [in] ntau - sequence number of scan point to be considered !! @param[out] dj_average - average of offset values in this run !! @param[out] dj_variance - varaiance of offset values in this run !! @param[out] dj_value - actual applied offset value !! !! @description !! This is user callable subroutine - public interface to fit_data module. subroutine Check_Dj_Offset(run,ntau, dj_average, dj_variance, dj_value) implicit none integer, intent(in) :: run integer, intent(in) :: ntau real(kind=DBL), intent(out) :: dj_average real(kind=DBL), intent(out) :: dj_variance real(kind=DBL), intent(out) :: dj_value integer :: i, nav dj_average = 0 dj_variance = 0 dj_value = 0 ! 040219 mm. nav = 0 do i=1, count_dj_offset if( runno_point_stat(1,i) == run ) then nav = nav + 1 dj_average = dj_average + dj_phaseoffsets(i) dj_variance = dj_variance + dj_phaseoffsets(i)**2 if( runno_point_stat(2,i) == ntau ) then dj_value = dj_phaseoffsets(i) endif endif enddo if( nav == 0) return dj_average = dj_average / nav dj_variance = sqrt(dj_variance/nav - dj_average**2) end subroutine Check_Dj_Offset end module fit_data sources/plot_utils.f90 +4 −4 Original line number Diff line number Diff line Loading @@ -338,7 +338,7 @@ contains ! choose between centerbin a if that seems not to be good, try to use central pixel Pixel => scan_point%centerbin if( Pixel%amplitude%value == 0d0 .or. Pixel%status .ne. PIXEL_OK) then !! if it is not present.. if( Pixel%amplitude%value == 0d0 .or. .not.is_pixel_ok(Pixel) ) then !! if it is not present.. Pixel => scan_point%pixelbin(ila,xpix,ypix) endif Loading Loading @@ -559,7 +559,7 @@ contains ! choose between centerbin a if that seems not to be good, try to use central pixel Pixel_ref => scan_point_ref%centerbin Pixel => scan_point%centerbin if( Pixel%amplitude%value == 0d0 .or. Pixel%status .ne. PIXEL_OK) then !! if it is not present.. if( Pixel%amplitude%value == 0d0 .or. .not.is_pixel_ok(Pixel) ) then !! if it is not present.. Pixel => scan_point%pixelbin(ila,xpix,ypix) Pixel_ref => scan_point_ref%pixelbin(ila,xpix,ypix) ! (PAZ) need to update the reference, bug #16 endif Loading Loading @@ -938,7 +938,7 @@ contains !! --> cross out invaid pixels if( Pixel%status .ne. PIXEL_OK) then if( .not.is_pixel_ok(Pixel) ) then ! call gr_setfillcolorind(106) ! call gr_fillarea(5,xhh,yhh) Loading Loading @@ -971,7 +971,7 @@ contains ipo: if(Pixel%status == PIXEL_OK) then ipo: if(is_pixel_ok(Pixel)) then call gr_setlinecolorind( 1 ) ! ..hope it is black else call gr_setlinecolorind( 2 ) Loading Loading
sources/fit_data00.f90deleted 100644 → 0 +0 −584 Original line number Diff line number Diff line !! ========================================================================== !! module fit_data !! ========================================================================== module fit_data use constants_module use data_types use echo_shapes use dump_data use fit_utils use mathutil use logger implicit none !!mm>> new to save offset info for reporting (experimental) integer, parameter :: MAX_TOTAL_POINTS = 10000 real(kind=DBL) :: DJOFFSET_LIMIT = 0.2d-6 integer :: count_dj_offset = 0 integer :: runno_point_stat(3,MAX_TOTAL_POINTS) = 0 real(kind=DBL) :: dj_phaseoffsets(MAX_TOTAL_POINTS) = 0 !!mm<< contains !! ======================================================================== !> @brief Fit echo data !! !! @param [in,out] data - scan_struct data (a @f$ \tau @f$ scan) !! @param [in] fit_flag - fitting level flag @see "phase_" flags in drspine_parameters !! @param [in,optional] phase_offset - initial phase offset !! !! @description !! This is user callable subroutine - public interface to fit_data module. subroutine fit_echo_data(data, fit_flag, phase_offset) type(scan_struct), target, intent(inout) :: data integer, intent(in) :: fit_flag real(kind=DBL), intent(in), optional :: phase_offset ! type(scan_data_struct), pointer :: sig_point, ref_point integer :: iscan real(kind=DBL) :: djoffset call msg_info('fit_echo_data', 'fitting echo data '//trim(data%file)//& ', role='//trim(cformat_role(data%role))//& ', flag='//trim(cformat_fit_flag(fit_flag))) djoffset = 0.0d0 if ( present(phase_offset) .and. data%number_of_points>0 ) & djoffset = get_field_integral(deg2rad(phase_offset), data%scan_point(1)%physics%lambda) !just need a wavelength ref_point => null() do iscan=1, data%number_of_points call msg_debug('fit_echo_data', msg_fmt("('scan point[0]: ', i4)", iscan)) sig_point => data%scan_point(iscan) if ( associated(sig_point%matching_ref) ) ref_point => sig_point%matching_ref !! call fit_center_bin(sig_point, ref_point, fit_flag) if ( iand(fit_flag, PHASE_FIT) /= 0 ) then ! fit phase for full b-w (no smoothing yet) call fit_echo_param(sig_point, ref_point, iand(fit_flag, not(PHASE_FIT_SMOOTH))) ! smooth if ( iand(fit_flag, PHASE_FIT_SMOOTH) /= 0 ) then call fit_echo_param(sig_point, sig_point, ior(PHASE_FIT, PHASE_FIT_SMOOTH)) endif if ( .not. associated(sig_point%matching_ref) ) sig_point%matching_ref => sig_point ref_point => sig_point endif if ( associated(sig_point%matching_ref) ) then call get_echo_param(sig_point, sig_point%matching_ref, .false.) ! extract if ( iand(fit_flag, PHASE_FIT_OFFSET) /=0 ) then if (iand(fit_flag, PHASE_FIT_RESET) /=0 ) djoffset=0 call get_phase_offset(sig_point, sig_point%matching_ref, djoffset) call get_echo_param(sig_point, sig_point%matching_ref, .true., djoffset) end if end if end do end subroutine fit_echo_data !! ======================================================================== !> @brief This subroutine fits phase scans (for resolution samples) !! !! @param dat_point - scan_data_struct pointer with phase scan data !! @param ref_point - scan_data_struct pointer with reference data flag (may be null) !! @param fit_flag - fitting level flag @see "phase_" flags in drspine_parameters !! !! @description !! This subroutine will perform symmetry phase scan fit over entire scan_data_struct !! (typically all detector pixels). !! It loops from the center of the detector in a spiral sequence. !! !! The supported flags are !! - PHASE_FIT - (implied) will perform phase fitting !! - PHASE_FIT_HOPPING - fit phase (full wavelength band) with "hopping" !! - PHASE_FIT_AMPPOS - fitted amplitude must be positive !! - PHASE_FIT_AMPNEG - fitted amplitude must be negative !! - PHASE_FIT_SMOOTH - smooth phase map and refit !! !! @note If ref_point is not null (i.e. associated) it will be used as a starting point !! for fitting (unless PHASE_FIT_SMOOTH is set). !! @note The ref_point may refer to the data itself. subroutine fit_echo_param(dat_point, ref_point, fit_flag) type(scan_data_struct),pointer :: dat_point, ref_point integer, intent(in) :: fit_flag ! !! it=0 fit phases only for the entire wavelength band integer, parameter :: it=0 ! fit full b-w integer :: ix, iy, ipix integer :: xpix, ypix type(lambda_struct) :: eshape_lams type(phase_scan_struct),pointer :: datpix, refpix real(kind=DBL) :: dj0 ! symmetry phase initial guess real(kind=DBL) :: ddj ! symmetry phase limits dj0 +/- ddj (typically +/- 90 degrees in phase) real(kind=DBL) :: cdj ! smoothing correction unit typically corresponding to 360 degree in phase ! integer, dimension(:), allocatable :: ixav integer, dimension(:), allocatable :: iyav !call print_image(97, scan_point%pixelbin(it,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAP1', z6)", fit_flag), .FALSE.) ddj = get_field_integral(reduction_parameters%symmetry_phase_max_dev, dat_point%physics%lambda) cdj = get_field_integral(reduction_parameters%symmetry_phase_corr_step, dat_point%physics%lambda) if ( iand(fit_flag, PHASE_FIT_SMOOTH)/=0 ) then call smooth_phase_map(dat_point%pixelbin(it,:,:)%delta_J_symm, & reduction_parameters%symmetry_phase_patch_size, ddj, cdj) !call print_image(97, dat_point%pixelbin(it,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAPS', z7)", fit_flag), .FALSE.) !return endif xpix = max(1, dat_point%no_xpix) ! at least one ypix = max(1, dat_point%no_ypix) allocate(ixav(xpix*ypix)) allocate(iyav(xpix*ypix)) call create_spiralad(ixav,iyav, xpix, ypix) call init_lambda_struct(eshape_lams, dat_point%spectrum(0)) pixel_loop: do ipix=1, xpix*ypix dj0 = 0.0_DBL refpix => null() !if (ipix==0) then ! datpix => dat_point%centerbin ! if ( associated(ref_point) ) refpix => ref_point%centerbin !else ix = ixav(ipix) iy = iyav(ipix) datpix => dat_point%pixelbin(it, ix, iy) if ( associated(ref_point) ) refpix => ref_point%pixelbin(it, ix, iy) !end if if ( iand(fit_flag, PHASE_FIT_SMOOTH)/=0 .and. ipix/=0 ) then dj0 = datpix%delta_J_symm%value else if ( associated(refpix) ) then dj0 = refpix%delta_J_symm%value endif call fit_phase_scan(datpix, eshape_lams, fit_flag, dj0, ddj) call msg_debug('fit_phase', clog_fit_params(datpix, 'pha:')) end do pixel_loop deallocate(ixav) deallocate(iyav) !call print_image(97, scan_point%pixelbin(it,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAP2', z6)", fit_flag), .FALSE.) end subroutine fit_echo_param !! ======================================================================== !> @brief This subroutine fits phase scans (for resolution samples) !! !! @param dat_point - scan_data_struct pointer with phase scan data !! @param ref_point - scan_data_struct pointer with reference data flag (may be null) !! @param fit_flag - fitting level flag @see "phase_" flags in drspine_parameters !! subroutine fit_center_bin(dat_point, ref_point, fit_flag) type(scan_data_struct),pointer :: dat_point, ref_point integer, intent(in) :: fit_flag ! type(lambda_struct) :: eshape_lams type(phase_scan_struct),pointer :: datpix, refpix real(kind=DBL) :: dj0 ! symmetry phase initial guess real(kind=DBL) :: ddj ! symmetry phase limits dj0 +/- ddj (typically +/- 90 degrees in phase) ! ddj = get_field_integral(reduction_parameters%symmetry_phase_max_dev, dat_point%physics%lambda) call init_lambda_struct(eshape_lams, dat_point%spectrum(0)) dj0 = 0.0_DBL refpix => null() datpix => dat_point%centerbin if ( associated(ref_point) ) then refpix => ref_point%centerbin dj0 = refpix%delta_J_symm%value end if call fit_phase_scan(datpix, eshape_lams, fit_flag, dj0, ddj) call msg_debug('fit_phase', clog_fit_params(datpix, 'pha:')) end subroutine fit_center_bin !! ============================================================================ !> @brief This subroutine gets (extracts) echo parameters using ref_point phase !! !! @param scan_point - scan_data_struct pointer with phase scan data !! @param ref_point - scan_data_struct pointer with reference data flag (cannot be null) !! !! @description !! !! @note The ref_point may refer to the data itself. subroutine get_echo_param(scan_point, res_point, apply_offset, offset) type(scan_data_struct) :: scan_point, res_point logical, intent(in) :: apply_offset real(kind=DBL), intent(in), optional :: offset ! integer :: it, ix, iy integer :: tpix, xpix, ypix type(lambda_struct) :: eshape_lams type(value_struct) :: dj real(kind=DBL) :: djoffset tpix = scan_point%no_lam xpix = scan_point%no_xpix ypix = scan_point%no_ypix call init_lambda_struct(eshape_lams, scan_point%spectrum(0)) djoffset = 0 if (apply_offset) then if (present(offset) ) djoffset = offset call msg_info('get_echo_param', & trim(msg_fmt("(1x,'applying phase offset: ',g12.5)", djoffset))//& trim(msg_fmt("(1x,'to scan[',i0,'/',i0,']:')", [scan_point%id, scan_point%point]))) !!mm>> new to save offset info for reporting (experimental) if(count_dj_offset < MAX_TOTAL_POINTS) count_dj_offset = count_dj_offset + 1 runno_point_stat(1,count_dj_offset) = scan_point%id runno_point_stat(2,count_dj_offset) = scan_point%point runno_point_stat(3,count_dj_offset) = int(abs(djoffset)/DJOFFSET_LIMIT) dj_phaseoffsets ( count_dj_offset) = djoffset !!mm<< end if do it=0, tpix !call msg_debug('get_echo_param', & ! msg_fmt("('scan point[',i2,']: ', i4)", [it ,iscan])) do ix=1, xpix do iy=1, ypix ! FIXME what if res and scan are one and the same? if ( res_point%pixelbin( 0, ix, iy)%status .ne. PIXEL_OK ) then scan_point%pixelbin(it, ix, iy)%status = & ior(scan_point%pixelbin(it, ix, iy)%status, & PIXEL_NO_RESOLUTION) cycle ! FIXME: shall we allow to fit without resolution?? endif dj = res_point%pixelbin(0, ix, iy)%delta_J_symm !if (apply_offset) then ! dj = add_value(dj, off) !end if !dj0 = res_point%pixelbin(0, ix, iy)%delta_J_symm%value !dje = sqrt(res_point%pixelbin(0, ix, iy)%delta_J_symm%sigma2) !call fit_phase_scan(scan_point%pixelbin(it, ix, iy), eshape_lams, PHASE_USE, dj0, dje) call fit_phase_scan(scan_point%pixelbin(it, ix, iy), eshape_lams, PHASE_USE, dj%value+djoffset, sqrt(dj%sigma2)) call msg_debug('get_echo_param', clog_fit_params(scan_point%pixelbin(it, ix, iy), 'ext:')) end do end do end do !call print_image(97, scan_point%pixelbin(0,:,:)%delta_J_symm , & ! msg_fmt("('PHASE_MAP', i4)", iscan), .FALSE.) end subroutine get_echo_param !! ======================================================================== !> @brief Perform actual symmetry phase fit !! !! @param phase_scan - phase scan data !! @param eshape_lams - wavelength spectrum !! @param fit_flag - fit level flag !! @param dj0 - initial guess for the symmetry phase !! @param djlim - phase minimization limits in Tm (Tesla-meters). !! @param tolerance - (optional) minimization tolerance !! !! @description !! Algorithm outline: !! - pre-fit check (only if fitting and not smoothing) !! - check if we have enough "wiggling" in the echo !! - check if the average in inside up/down bounds !! - fit the phase withing the limits (if PHASE_FIT is set) !! - extract echo parameters (average, amplitude, phase) with errors !! - post-fit check !! - check @f$\chi^2@f$ !! - check fit amplitude and average !! subroutine fit_phase_scan(phase_scan, eshape_lams, fit_flag, dj0, dj0_lim, tolerance) type(phase_scan_struct), intent(inout) :: phase_scan ! input data type(lambda_struct), intent(inout) :: eshape_lams ! wavelength spectrum integer, intent(in) :: fit_flag ! fit flag real(kind=DBL), intent(in) :: dj0 ! initial guess real(kind=DBL), intent(in) :: dj0_lim ! fit limits real(kind=DBL), intent(in), optional :: tolerance ! fit tolerance ! integer, parameter :: NHOPS = 3 real(kind=DBL), dimension(-NHOPS:NHOPS,2) :: dj ! dj bracket(s) real(kind=DBL) :: dj0_fit, dj0_err ! fitted phase and error real(kind=DBL) :: dj0_tol ! tolerance real(kind=DBL) :: chi2_fit ! absolute chi2 minimum real(kind=DBL) :: xchi, xdj ! tmp vars integer :: stat ! minimization status integer :: i real(kind=DBL), dimension(2) :: dj0_hard ! type(value_struct) :: amp ! optional parameters dj0_tol = reduction_parameters%symmetry_phase_tolerance if ( present(tolerance) ) dj0_tol = tolerance if(phase_scan%amplitude%state == VAL_UNDEFINED) then ! (PAZ: FIXME, weird way to test for "not checked" raw data) if ( iand(fit_flag, PHASE_FIT) > 0 .and. iand(fit_flag, PHASE_FIT_SMOOTH)==0 ) then call check_raw_data(phase_scan, reduction_parameters%min_sigma_relative) end if endif if (phase_scan%status/=PIXEL_OK) return eshape_lams%nbin = [phase_scan%tbin_1, phase_scan%tbin_2] dj0_fit = dj0 dj0_err = HUGE(dj0) chi2_fit = HUGE(dj0) dj0_hard(1) = MINVAL(phase_scan%delta_J(1:phase_scan%n_down)%value) dj0_hard(2) = MAXVAL(phase_scan%delta_J(1:phase_scan%n_down)%value) if ( iand(fit_flag, PHASE_FIT) /= 0 ) then if ( iand(fit_flag, PHASE_FIT_HOPPING) == 0 ) then dj(0,1) = MAX(dj0_hard(1), dj0 - dj0_lim) ! brent needs a bracket dj(0,2) = MIN(dj0_hard(2), dj0 + dj0_lim) if ( dj(0,1)>dj(0,2) ) call swap(dj(0,1), dj(0,2)) stat = linextpha1(eshape, phase_scan, eshape_lams, dj(0,:), dj0_tol) !minimize if (stat<0) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED) return end if dj0_fit = 0.5*(dj(0,1)+dj(0,2)) else do i=-NHOPS,NHOPS dj(i,1) = MAX(dj0_hard(1), 2*(i-1)*dj0_lim) ! brent needs a bracket dj(i,2) = MIN(dj0_hard(2), 2*(i+1)*dj0_lim) end do ! fit with 'hopping' stat = linextpha2(eshape, phase_scan, eshape_lams, dj, dj0_tol) !minimize if (stat<0) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED) return end if hops: do i=-NHOPS,NHOPS xdj = 0.5*(dj(i,1)+dj(i,2)) call linextract(eshape, phase_scan, eshape_lams, xdj, xchi, amp=amp) if (iand(fit_flag, PHASE_FIT_AMPPOS)/=0 .and. amp%value<0) cycle hops if (iand(fit_flag, PHASE_FIT_AMPNEG)/=0 .and. amp%value>0) cycle hops if ( xchi<chi2_fit ) then dj0_fit = xdj chi2_fit = xchi end if end do hops end if end if if (phase_scan%status/=PIXEL_OK) return call get_phase_fit_values(eshape, phase_scan, eshape_lams, dj0_fit) ! get fit values phase_scan%fqt = get_fqt(phase_scan) call check_fit_par(phase_scan, fit_flag, & reduction_parameters%fqt_min_value,& reduction_parameters%average_max_deviation, & reduction_parameters%max_chisquare) end subroutine fit_phase_scan ! ======================================================================== !> Helper subroutine. !! Check if we have enough "wiggling" in the echo signal by checking !! simga/average of the data against a threshold !! Use this function *before* fitting phase !! @param phase_scan - phase_scan_struct data !! @param threshold - min value for sigma/average subroutine check_raw_data(phase_scan, threshold) type(phase_scan_struct), intent(inout) :: phase_scan real(kind=DBL), intent(in) :: threshold ! real(kind=DBL) :: ave, sig, up, dn ! ave = phase_scan%average%value sig = sqrt(phase_scan%average%sigma2) !! up = max(phase_scan%up%value, phase_scan%down%value) !! dn = min(phase_scan%up%value, phase_scan%down%value) up = phase_scan%up%value dn = phase_scan%up%value ! check counting statistics if ( phase_scan%average_raw%value.lt.reduction_parameters%min_counts_per_pixel) then phase_scan%status = ior(phase_scan%status, PIXEL_STATISTICS) !call msg_debug('check_raw_data', trim(phase_scan%cpixel)//' not enough statistics'//& ! trim(msg_fmt('(" counts=",g12.6)', phase_scan%average_raw%value))//& ! trim(msg_fmt('(" status=",i0 )', phase_scan%status ))) return end if ! if (ave<0) then phase_scan%status = ior(phase_scan%status, PIXEL_AMPLITUDE_RAW) else if (sig/ave<threshold) & phase_scan%status = ior(phase_scan%status, PIXEL_AMPLITUDE_RAW) if (ave<dn .or. up<ave) & phase_scan%status = ior(phase_scan%status, PIXEL_AVERAGE_OUT_OF_BOUNDS) end if end subroutine check_raw_data ! ======================================================================== !> Helper subroutine. !! 1. Check if the relative amplitude of the fit (udratio) is high enough !! 2. Check if the fit average is consistent with up/down average !! !! Use this function *after* fitting phase !! @param phase_scan - phase_scan_struct data !! @param fit_flag - fitting flag !! @param threshold_amp - min value for sigma/average !! @param threshold_ave - max relative deviation !! @param threshold_chi2 - max allowed chi2 subroutine check_fit_par(phase_scan, fit_flag, threshold_amp, threshold_ave, threshold_chi2) type(phase_scan_struct), intent(inout) :: phase_scan integer, intent(in) :: fit_flag real(kind=DBL), intent(in) :: threshold_amp real(kind=DBL), intent(in) :: threshold_ave real(kind=DBL), intent(in) :: threshold_chi2 ! real(kind=DBL) :: ave, amp, fqt, uda ave = phase_scan%average%value amp = phase_scan%amplitude%value fqt = phase_scan%fqt%value uda = (phase_scan%up%value + phase_scan%down%value)/2 if ( phase_scan%chi_squared > threshold_chi2 ) & phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_CHISQUARED) if ( iand(fit_flag, PHASE_FIT_AMPPOS) > 0 .and. amp<0 ) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_NEGATIVE) endif if ( iand(fit_flag, PHASE_FIT_AMPNEG) > 0 .and. amp>0 ) then phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_POSITIVE) endif if ( iand(fit_flag, PHASE_FIT) > 0 ) then if ( ave>0 .and. abs(fqt)<threshold_amp) & phase_scan%status = ior(phase_scan%status, PIXEL_FIT_FAILED_AMPLITUDE) if ( abs(ave-uda)/uda > threshold_ave ) & phase_scan%status = ior(phase_scan%status, PIXEL_AVERAGE_INCONSISTENT) endif return end subroutine check_fit_par ! ======================================================================== !> Helper subroutine. !! Log fit parameters function clog_fit_params(phase_scan, cprefix) result(cline) character(len=MAX_LINE_LENGTH) :: cline type(phase_scan_struct), intent(in) :: phase_scan character(len=*), intent(in) :: cprefix write(cline,'(a,1x,a,2(1x,g13.6),4(1x,a),1x,g13.6,z8)') & trim(cprefix), & trim(phase_scan%cpixel), & !indices(1), indices(2), & phase_scan%Q_abs*ANGSTROEM , & phase_scan%tau/NS, & trim(cformat_value(phase_scan%average, 'g13.6', MAX_NAME_LENGTH)),& trim(cformat_value(phase_scan%amplitude,'g13.6', MAX_NAME_LENGTH)),& trim(cformat_value(phase_scan%delta_J_symm, 'g13.6', MAX_NAME_LENGTH)),& trim(cformat_value(phase_scan%fqt, 'g13.6', MAX_NAME_LENGTH)),& phase_scan%chi_squared, & phase_scan%status return end function clog_fit_params function select_best_phase(sfun, phase_scan, eshape_lams, dJ0) result(res) procedure(echo_shape_function) :: sfun type(phase_scan_struct), intent(in) :: phase_scan type(lambda_struct), intent(in) :: eshape_lams real(kind=DBL), intent(in), dimension(:) :: dJ0 real(kind=DBL) :: res ! real(kind=DBL) :: chi2, chi2_min type(value_struct) :: amp integer :: k ! res = 0.0_DBL chi2_min = HUGE(1.0_DBL) !print *, phase_scan%delta_J(1)%value, phase_scan%delta_J(phase_scan%n_down)%value do k=1, SIZE(dJ0) call linextract(sfun, phase_scan, eshape_lams, dJ0(k), chi2, amp=amp) !print *, k, dJ0(k), chi2, amp%value if ( dJ0(k) < phase_scan%delta_J(1)%value .or. phase_scan%delta_J(phase_scan%n_down)%value < dJ0(k)) cycle if ( amp%value<0) cycle !print *, 'OK' if ( chi2 < chi2_min ) then res = dJ0(k) chi2_min = chi2 end if end do end function select_best_phase !!mm>> !! ======================================================================== !!> @brief Check_Dj_Offset(run,ntau, dj_variance, dj_value) !! !! @param [in] run - numor of run to be considered !! @param [in] ntau - sequence number of scan point to be considered !! @param[out] dj_average - average of offset values in this run !! @param[out] dj_variance - varaiance of offset values in this run !! @param[out] dj_value - actual applied offset value !! !! @description !! This is user callable subroutine - public interface to fit_data module. subroutine Check_Dj_Offset(run,ntau, dj_average, dj_variance, dj_value) implicit none integer, intent(in) :: run integer, intent(in) :: ntau real(kind=DBL), intent(out) :: dj_average real(kind=DBL), intent(out) :: dj_variance real(kind=DBL), intent(out) :: dj_value integer :: i, nav dj_average = 0 dj_variance = 0 dj_value = 0 ! 040219 mm. nav = 0 do i=1, count_dj_offset if( runno_point_stat(1,i) == run ) then nav = nav + 1 dj_average = dj_average + dj_phaseoffsets(i) dj_variance = dj_variance + dj_phaseoffsets(i)**2 if( runno_point_stat(2,i) == ntau ) then dj_value = dj_phaseoffsets(i) endif endif enddo if( nav == 0) return dj_average = dj_average / nav dj_variance = sqrt(dj_variance/nav - dj_average**2) end subroutine Check_Dj_Offset end module fit_data
sources/plot_utils.f90 +4 −4 Original line number Diff line number Diff line Loading @@ -338,7 +338,7 @@ contains ! choose between centerbin a if that seems not to be good, try to use central pixel Pixel => scan_point%centerbin if( Pixel%amplitude%value == 0d0 .or. Pixel%status .ne. PIXEL_OK) then !! if it is not present.. if( Pixel%amplitude%value == 0d0 .or. .not.is_pixel_ok(Pixel) ) then !! if it is not present.. Pixel => scan_point%pixelbin(ila,xpix,ypix) endif Loading Loading @@ -559,7 +559,7 @@ contains ! choose between centerbin a if that seems not to be good, try to use central pixel Pixel_ref => scan_point_ref%centerbin Pixel => scan_point%centerbin if( Pixel%amplitude%value == 0d0 .or. Pixel%status .ne. PIXEL_OK) then !! if it is not present.. if( Pixel%amplitude%value == 0d0 .or. .not.is_pixel_ok(Pixel) ) then !! if it is not present.. Pixel => scan_point%pixelbin(ila,xpix,ypix) Pixel_ref => scan_point_ref%pixelbin(ila,xpix,ypix) ! (PAZ) need to update the reference, bug #16 endif Loading Loading @@ -938,7 +938,7 @@ contains !! --> cross out invaid pixels if( Pixel%status .ne. PIXEL_OK) then if( .not.is_pixel_ok(Pixel) ) then ! call gr_setfillcolorind(106) ! call gr_fillarea(5,xhh,yhh) Loading Loading @@ -971,7 +971,7 @@ contains ipo: if(Pixel%status == PIXEL_OK) then ipo: if(is_pixel_ok(Pixel)) then call gr_setlinecolorind( 1 ) ! ..hope it is black else call gr_setlinecolorind( 2 ) Loading
