Commit 2d819bf3 authored by Michael Monkenbusch(AB)'s avatar Michael Monkenbusch(AB)
Browse files

averaging of matching refs and bgrs in reduce_data calculate_sqt_raw, still... 

averaging of matching refs and bgrs in reduce_data calculate_sqt_raw, still contains test writes flagged with !!TEST!!
parent ed5398ce

sources/base_types.f90

+24 −0
Original line number Diff line number Diff line
@@ -95,6 +95,7 @@ module base_types 


  !> weighted mean

  interface ave_value

     module procedure ave_value0    !>>new>>:mm

     module procedure ave_value1

     module procedure ave_value2

     module procedure ave_value3
@@ -424,6 +425,29 @@ contains 
    res%value_struct = div_value2(arg1%value_struct, arg2)

  end function div_namedvalue2



  !> combine two values with error weights

  function ave_value0(arg1, arg2) result(res)

    type(value_struct), intent(in) :: arg1

    type(value_struct), intent(in) :: arg2

    type(value_struct) :: res

    real(kind=DBL)     :: weight1=0.0_DBL, weight2=0.0_DBL



    if(arg1%sigma2>0.0_DBL) weight1 = 1/arg1%sigma2

    if(arg2%sigma2>0.0_DBL) weight2 = 1/arg2%sigma2



    if( weight1+weight2 > 0.0_DBL)  then 

       res%state  = VAL_WITH_ERROR

    else

       call init_value(res)

       return

    endif



    res%value  = (arg1%value * weight1 + arg2%value * weight2)/(weight1+weight2)  

    res%sigma2 = 1.0_DBL / (weight1+weight2)

!!TEST!! write(*,'(a,3(i3,2f18.9))')"TAVE: ", arg1, arg2, res 

!!TEST!! write(*,'(4(a,e15.7))')"TAVE w1=",weight1," w2=",weight2," a1=",arg1%value," a2=",arg2%value

  

  end function ave_value0





 !> calculate weighted average

sources/drspine_parameters.f90

+1 −0
Original line number Diff line number Diff line
@@ -86,6 +86,7 @@ module drspine_parameters 
  integer, parameter      :: PIXEL_NO_RESOLUTION         = INT(Z'1000') ! no resolution for pixel

  integer, parameter      :: PIXEL_NO_BACKGROUND         = INT(Z'2000') ! no background for pixel

  integer, parameter      :: PIXEL_OFFSET_FAILED         = INT(Z'4000') ! offset optimization failed

  integer, parameter      :: PIXEL_COMBINATION_FAILED    = INT(Z'8000') ! offset optimization failed



  !

  integer, parameter      :: PHASE_USE           = INT(Z'000000')  ! use given phase (do not fit)

sources/matching.f90

+10 −10
Original line number Diff line number Diff line
@@ -192,7 +192,7 @@ contains 
          if(.not. associated( source_data%scan_point(i)%matching_bgr)) &

              call clear_matches(source_data%scan_point(i)%matching_bgr_arr)

      endif

write(*,*)"TESTTAUMATCH 1 : clear  ", i   !!TEST!!

!!TEST!! write(*,*)"TESTTAUMATCH 1 : clear  ", i   !!TEST!!

 ! <<new <<

       target_loop: do j=1, target_data%number_of_points

          if (taupoint_match(source_data%scan_point(i), target_data%scan_point(j))) then
@@ -204,15 +204,15 @@ write(*,*)"TESTTAUMATCH 1 : clear ", i !!TEST!! 
                                      source_data%scan_point(i)%matching_ref  ) 

!!TEST!!                 write(*,*)"TESTTAUMATCH 2 : ref  ", i, j,  target_data%scan_point(j)%id, &   !!TEST!!

!!TEST!!                           source_data%scan_point(i)%matching_ref_arr(1)%ptr%id               !!TEST!!

!>TEST

                do l=1,size(source_data%scan_point(i)%matching_ref_arr)

                  if(.not.associated(source_data%scan_point(i)%matching_ref_arr(l)%ptr)) then

                    write(*,*) l, " not assoc"

                  else

                    write(*,*) l,source_data%scan_point(i)%matching_ref_arr(l)%ptr%id 

                  endif

                enddo

!>TEST

!!!TEST!!>TEST

!!TEST!!        do l=1,size(source_data%scan_point(i)%matching_ref_arr)

!!TEST!!          if(.not.associated(source_data%scan_point(i)%matching_ref_arr(l)%ptr)) then

!!TEST!!            write(*,*) l, " not assoc"

!!TEST!!          else

!!TEST!!            write(*,*) l,source_data%scan_point(i)%matching_ref_arr(l)%ptr%id 

!!TEST!!          endif

!!TEST!!        enddo

!!!TEST!!>TEST

 ! <<new<<

                cmsg = 'res'

             else if(has_role(target_data%role, ROLE_BACKGROUND)) then

sources/reduce_data.f90

+100 −13
Original line number Diff line number Diff line
@@ -83,12 +83,14 @@ contains 
    !

    integer :: iscan

    type(value_struct)      :: fqt

    type(phase_scan_struct) :: resol

    type(phase_scan_struct) :: bgr

    integer :: it, ix, iy



    logical :: has_bgr

    real(kind=DBL) :: tfac



    do iscan=1, sig_data%number_of_points

ds: do iscan=1, sig_data%number_of_points

       has_bgr = .false.

       sig: associate(sig_point => sig_data%scan_point(iscan))

       if ( .not. associated(sig_point%matching_ref) ) cycle
@@ -100,14 +102,23 @@ contains 
       end if

       res: associate(res_point => sig_point%matching_ref,&

                      bgr_point => sig_point%matching_bgr)

       do it=0, sig_point%no_lam

          do ix=1, sig_point%no_xpix

             do iy=1, sig_point%no_ypix

dit:   do it=0, sig_point%no_lam

dix:      do ix=1, sig_point%no_xpix

diy:         do iy=1, sig_point%no_ypix

                if ( sig_point%pixelbin(it, ix, iy)%status /= PIXEL_OK ) cycle



if(ix==4 .and. iy==4) write(*,*) &  !!TEST!!

"T combine ref",res_point%id,sig_point%matching_ref_arr(1)%ptr%id,& !!TEST!!

 associated(sig_point%matching_ref, sig_point%matching_ref_arr(1)%ptr)



call combine_runs(sig_point%matching_ref_arr, it, ix, iy, resol  ) ! later use resol%fqt

                if ( res_point%pixelbin(it, ix, iy)%status /= PIXEL_OK ) cycle

                if ( bgr_flag > 0 .and. has_bgr) then

if(ix==4 .and. iy==4) write(*,*)"T combine bgr" !!TEST!!

 call combine_runs(sig_point%matching_bgr_arr, it, ix, iy, bgr )    ! for later use instead of bgr_point%pixelbin

                   if (bgr_point%pixelbin(it, ix, iy)%status /= PIXEL_OK ) then

                      sig_point%pixelbin(it, ix, iy)%status = ior(sig_point%pixelbin(it, ix, iy)%status, PIXEL_NO_BACKGROUND)

                      sig_point%pixelbin(it, ix, iy)%status = & 

                          ior(sig_point%pixelbin(it, ix, iy)%status, PIXEL_NO_BACKGROUND)

                      cycle

                   end if

                   if ( transmission_factor > 0 ) then
@@ -115,23 +126,99 @@ contains 
                   else

                        tfac = sig_point%pixelbin(it, ix, iy)%transmission_factor%value

                   end if

                   fqt = get_fqt_bgr(sig_point%pixelbin(it, ix, iy), bgr_point%pixelbin(it, ix, iy), tfac, volume_fraction)



!!TEST!! 

                   fqt = get_fqt_bgr(sig_point%pixelbin(it, ix, iy), bgr_point%pixelbin(it, ix, iy), tfac, volume_fraction) !!TEST!! 

if(ix==4 .and. iy==4) write(*,*)"T fqt1_sngl(old)",fqt !!TEST!!





                   fqt = get_fqt_bgr(sig_point%pixelbin(it, ix, iy), bgr, tfac, volume_fraction)

if(ix==4 .and. iy==4) write(*,*)"T fqt2_mult(new)",fqt !!TEST!!

                else

                   fqt = sig_point%pixelbin(it, ix, iy)%fqt

                endif

                sig_point%pixelbin(it, ix, iy)%fqtRes     = div_value(fqt, res_point%pixelbin(it, ix, iy)%fqt)

                sig_point%pixelbin(it, ix, iy)%resolution = res_point%pixelbin(it, ix, iy)%fqt

              end do

           end do

        end do

 



!!TEST!! if(ix==4 .and. iy==4) write(*,*) & !!TEST!!

!!TEST!! "T resA",sig_point%matching_ref%id,sig_point%matching_ref%point, & !!TEST!!

!!TEST!! sig_point%matching_ref%pixelbin(it, ix, iy)%fqt, resol%fqt !!TEST!!

!!TEST!! if(ix==4 .and. iy==4) write(*,*) & !!TEST!!

!!TEST!! "T resB",res_point%id,res_point%point, & !!TEST!!

!!TEST!! res_point%pixelbin(it, ix, iy)%fqt, resol%fqt !!TEST!!

!? resol%fqt = res_point%pixelbin(it, ix, iy)%fqt



!?                sig_point%pixelbin(it, ix, iy)%fqtRes     = div_value(fqt, res_point%pixelbin(it, ix, iy)%fqt)

!?                sig_point%pixelbin(it, ix, iy)%resolution = res_point%pixelbin(it, ix, iy)%fqt

                sig_point%pixelbin(it, ix, iy)%fqtRes     = div_value(fqt, resol%fqt)

                sig_point%pixelbin(it, ix, iy)%resolution = resol%fqt

              end do diy

           end do dix

        end do dit

        end associate res

        end associate sig

    end do

    end do ds



    !close(iunit)



  end subroutine calculate_sqt_raw



  !>@brief

  !! combine equivalent pixels from different runs

  subroutine combine_runs(run_pointer_array, it, ix, iy, averages )    

     type(scan_data_struct_pointer), dimension(:), intent(in) :: run_pointer_array

     integer                                     , intent(in) :: it, ix, iy        ! pixel addresses

     type(phase_scan_struct)                     , intent(out):: averages 



     integer  :: i, stat

     type(value_struct) :: amp, ave, up, down, dj, fqt



     call init_phase_scan_struct(averages, 1)    !! init necessary ??

 

     stat =  PIXEL_COMBINATION_FAILED 



     do i=1, size(run_pointer_array)

       if(.not. associated( run_pointer_array(i)%ptr))  cycle  !! or more restrictive but faster exit instead of cycle



!if(ix==4 .and. iy==4) write(*,*)"Tcombine_runs: ",it,ix,iy," pt: ",i,run_pointer_array(i)%ptr%id  !!TEST!! 



       if( run_pointer_array(i)%ptr%pixelbin(it,ix,iy)%status /= PIXEL_OK ) cycle  

       if(i == 1 ) then

         amp   =  run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%amplitude  

         ave   =  run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%average 

         up    =  run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%up 

         down  =  run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%down 

         dj    =  run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%delta_j_symm 

         fqt   =  run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%fqt                

       else

         amp   = ave_value( amp , run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%amplitude ) 

         ave   = ave_value( ave , run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%average )

         up    = ave_value( up  , run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%up )

         down  = ave_value( down, run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%down )

         dj    = ave_value( dj  , run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%delta_j_symm )

         fqt   = ave_value( fqt , run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%fqt  )           !! ??

       endif

       stat  = PIXEL_OK



if(ix==4 .and. iy==4) write(*,'(a,i4,i9,2(i3,2f18.9))') "Tcombine_runs i amp fqt: ",i,run_pointer_array(i)%ptr%id, amp,fqt !!TEST!!

!!TEST!! if(ix==4 .and. iy==4) write(*,'(a,i4,i9,2(i3,2f18.9))') "Tcombine_runs o amp fqt: ",i,run_pointer_array(i)%ptr%id,&!!TEST!! 

!!TEST!! run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%amplitude ,  run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%fqt !!TEST!! 

!!TEST!! if(ix==4 .and. iy==4) write(*,*)"FQT ", 2*amp%value/(up%value-down%value), run_pointer_array(i)%ptr%pixelbin(it, ix, iy)%fqt !!TEST!! 

     enddo



     averages%amplitude            =  amp  

     averages%average              =  ave  

     averages%up                   =  up    

     averages%down                 =  down 

     averages%delta_j_symm         =  dj   

     averages%fqt                  =  fqt       

     averages%status               =  stat       

    

 !!TEST!!  if(ix==4 .and. iy==4) write(*,'(a,i4,i9,2(i3,2f18.9))')"Tcombine_runs R amp fqt: ",i,run_pointer_array(1)%ptr%id,amp,fqt !!TEST!!



 

  end subroutine combine_runs

 



  !>@brief

  !! write collection header

  subroutine collect_sqt_raw_header(iunit)