implemented xrpd fit

VERSION

@@ -1 +1 @@
-1.2.312
+1.2.313

src/crystalmod@powderproc.F90

@@ -956,7 +956,7 @@ module subroutine xrpd_peaks_from_profile_file(p,xyfile,rms,errmsg,verbose0,&
 
     ! bail out if NLOPT is not available
 #ifndef HAVE_NLOPT
-    errmsg = "gaussian_compare can only be used if nlopt is available"
+    errmsg = "xrpd_peaks_from_profile_file can only be used if nlopt is available"
     rms = huge(1d0)
     return
 #else
@@ -965,7 +965,7 @@ module subroutine xrpd_peaks_from_profile_file(p,xyfile,rms,errmsg,verbose0,&
     integer, allocatable :: pid(:), io(:)
     real*8, allocatable :: x(:), y(:), pth2(:), phei(:), prm(:), lb(:), ub(:), yfit(:), ysum(:)
     real*8, allocatable :: yread(:)
-    real*8 :: fac, minx, maxx, maxy, ssq, maxa, xdif, xshift
+    real*8 :: fac, minx, maxx, maxy, ssq, maxa, xdif, xshift, ymean
     logical :: ok
     integer :: npeaks, npeaks_
     integer*8 :: opt
@@ -1008,6 +1008,8 @@ module subroutine xrpd_peaks_from_profile_file(p,xyfile,rms,errmsg,verbose0,&
 
     ! if default ymax_detect, use the median of the data
     if (ymax_detect == -huge(1d0)) then
+       ymean = sum(y - minval(y))/size(y,1)
+
        ! straightforward calculation of the median
        ysum = y
        allocate(io(n))

src/struct_drivers@proc.F90

@@ -4842,15 +4842,18 @@ end subroutine struct_bz
   !> Tools for editing and processing 1D x-ray powder diffraction
   !> patterns.
   module subroutine struct_xrpd(line0)
-    use crystalmod, only: david_sivia_calculate_background
+    use crystalmod, only: david_sivia_calculate_background, xrpd_peaklist
+    use global, only: fileroot
     use tools_io, only: uout, getword, lgetword, equal, ferror, faterr, isinteger, string,&
-       file_read_xy, fopen_write, fclose
+       file_read_xy, fopen_write, fclose, isreal
     character*(*), intent(in) :: line0
 
     character(len=:), allocatable :: word, xyfile, file, errmsg
-    integer :: lp, nknot, n, lu, i
+    integer :: lp, nknot, n, lu, i, nadj
     logical :: ok
-    real*8, allocatable :: x(:), y(:), yout(:)
+    real*8, allocatable :: x(:), y(:), yout(:), yfit(:)
+    real*8 :: ymax_peakdetect, rms
+    type(xrpd_peaklist) :: p
 
     integer, parameter :: nknot_def = 20
 
@@ -4886,15 +4889,50 @@ module subroutine struct_xrpd(line0)
 
        ! write the background to the file
        lu = fopen_write(file)
-       write (lu,'("## x  ybackground  yobs-ybackground  yobs")')
+       write (lu,'("## x  yobs-ybackground  ybackground  yobs")')
        do i = 1, n
           write (lu,'(4(A,X))') string(x(i),'f',decimal=10),&
-             string(yout(i),'e',decimal=10), string(y(i)-yout(i),'e',decimal=10),&
+             string(y(i)-yout(i),'e',decimal=10), string(yout(i),'e',decimal=10),&
              string(y(i),'e',decimal=10)
        end do
        call fclose(lu)
+
     elseif (equal(word,"fit")) then
-       write (*,*) "fixme2"
+       ! fit an experimental XRPD profile
+       !   FIT file-xy.s [ymax_peakdetect.r] [nadj.i]
+
+       ! header
+       write (uout,'("+ FIT: Fitting a peak profile to experimental XRPD data")')
+
+       ! read the input options
+       file = getword(line0,lp)
+       ymax_peakdetect = -huge(1d0)
+       ok = isreal(ymax_peakdetect,line0,lp)
+       if (ok) then
+          nadj = 2
+          ok = isinteger(nadj,line0,lp)
+          if (nadj /= 1 .and. nadj /= 2) &
+             call ferror('struct_xrpd','Number of adjacent points (NADJ) must be 1 or 2',faterr)
+       end if
+
+       ! run the peak fit
+       call p%from_profile_file(file,rms,errmsg,.true.,ymax_peakdetect,nadj,xorig=x,yorig=y,ycalc=yfit)
+       if (len_trim(errmsg) > 0) &
+          call ferror('struct_xrpd',errmsg,faterr)
+
+       ! write the profile to disk
+       lu = fopen_write(fileroot // "_fit.dat")
+       write (lu,'("## x yorig ycalc")')
+       do i = 1, size(x,1)
+          write (lu,'(3(A," "))') string(x(i),'f',decimal=10), string(y(i),'f',decimal=10),&
+             string(yfit(i),'f',decimal=10)
+       end do
+       call fclose(lu)
+
+       ! final list of peaks to disk
+       file = fileroot // ".peaks"
+       write (uout,'("+ List of peaks written to file: ",A)') file
+       call p%write(file)
     elseif (equal(word,"refit")) then
        write (*,*) "fixme3"
     else

src/tricks.F90

@@ -32,7 +32,6 @@ module tricks
   private :: trick_bfgs
   private :: trick_compare_deformed
   private :: trick_check_valence
-  private :: trick_profile_fit
   private :: trick_profile_refit
 
 contains
@@ -59,8 +58,6 @@ subroutine trick(line0)
        call trick_compare_deformed(line0(lp:))
     else if (equal(word,'check_valence')) then
        call trick_check_valence(line0(lp:))
-    else if (equal(word,'profile_fit')) then
-       call trick_profile_fit(line0(lp:))
     else if (equal(word,'profile_refit')) then
        call trick_profile_refit(line0(lp:))
     else
@@ -2914,75 +2911,6 @@ subroutine trick_check_valence(line0)
 
   end subroutine trick_check_valence
 
-  ! PROFILE_FIT file-xy.s ymax_peakdetect.r nadj.i
-  subroutine trick_profile_fit(line0)
-#ifndef HAVE_NLOPT
-    use tools_io, only: ferror, faterr, uout, fopen_write
-    character*(*), intent(in) :: line0
-
-    call ferror("trick_profile_fit","trick_profile_fit can only be used if nlopt is available",faterr)
-#else
-    use crystalmod, only: xrpd_peaklist
-    use tools_io, only: ferror, faterr, uout, getword, isreal, isinteger, string, fopen_write, fclose,&
-       string
-    character*(*), intent(in) :: line0
-
-    integer :: lpo
-    character(len=:), allocatable :: file, errmsg
-    real*8, allocatable :: x(:), y(:), yfit(:)
-    real*8 :: ymax_peakdetect, rms
-    integer :: nadj, lu, i, n
-    logical :: ok
-    type(xrpd_peaklist) :: p
-
-    ! header
-    write (uout,'("* Trick: profile fit")')
-
-    ! read the input options
-    lpo = 1
-    file = getword(line0,lpo)
-    ymax_peakdetect = -huge(1d0)
-    ok = isreal(ymax_peakdetect,line0,lpo)
-    if (ok) then
-       nadj = 2
-       ok = isinteger(nadj,line0,lpo)
-       if (nadj /= 1 .and. nadj /= 2) &
-          call ferror('trick','Number of adjacent points must be 1 or 2',faterr)
-    end if
-
-    ! run the peak fit
-    call p%from_profile_file(file,rms,errmsg,.true.,ymax_peakdetect,nadj,xorig=x,yorig=y,ycalc=yfit)
-    if (len_trim(errmsg) > 0) &
-       call ferror('trick_profile_fit',errmsg,faterr)
-
-    ! write the profile to disk
-    lu = fopen_write("fit.dat")
-    write (lu,'("## x yorig ycalc")')
-    do i = 1, size(x,1)
-       write (lu,'(3(A," "))') string(x(i),'f',decimal=10), string(y(i),'f',decimal=10),&
-          string(yfit(i),'f',decimal=10)
-    end do
-    call fclose(lu)
-
-    ! calculate final profile and write it to disk
-    n = 3000
-    call p%calculate_profile(n,x,y,errmsg)
-    if (len_trim(errmsg) > 0) &
-       call ferror('trick_profile_fit',errmsg,faterr)
-    lu = fopen_write("fitext.dat")
-    write (lu,'("## x ycalc")')
-    do i = 1, n
-       write (lu,'(2(A," "))') string(x(i),'f',decimal=10), string(y(i),'f',decimal=10)
-    end do
-    call fclose(lu)
-
-    ! final list of peaks to disk
-    write (uout,'("+ List of peaks written to file: fit.peaks")')
-    call p%write("fit.peaks")
-
-#endif
-  end subroutine trick_profile_fit
-
   ! PROFILE_REFIT file-xy.s fit.peaks
   subroutine trick_profile_refit(line0)
     use crystalmod, only: xrpd_peaklist
@@ -3037,7 +2965,7 @@ subroutine trick_profile_refit(line0)
     n = 3000
     call p%calculate_profile(n,x,y,errmsg)
     if (len_trim(errmsg) > 0) &
-       call ferror('trick_profile_fit',errmsg,faterr)
+       call ferror('trick_profile_refit',errmsg,faterr)
     lu = fopen_write("fitext.dat")
     write (lu,'("## x ycalc")')
     do i = 1, n