program TriCS_CCL c ----------------- c c fits data from ccl-files and write *.col file for upals c and rafin.dat file for position refining c c replaces "none" c 8.5.99/21.6.99/15.7.99 Juerg Schefer c 29.11.99: Lorentz-Korrektur SJ33 c 10.4.2000: flat background SJ33 c 31.10.2000: transfered to UNIX ZO33 c 22.11.2000: problem with printing is fixed c 02.02.01: h,k,l may be real ZM33 c 17.06.2010: Lorentz correction for tilt Geometry, SJ33 c c link ccl,fit4_shr/opt ! attention: check for last fit_n c implicit none integer np_max parameter (np_max=1000) real output(20),th2,omega_cal,omega_fit,chi,phi,yint,yy1, * dint,delta,yoverall,h,k,l,s,temperature,range real omega_err real output2(20), xx(np_max),lor_cor,fwhm,radian,degrees real bgr, wmult logical plot,iexp,flat_backgr integer lorentz,comens ! SJ33, June 2010 integer next_flag integer iall,iincom,icom,icenter character*1 answer,hard character*128 input_files character*128 output_file1 character*128 output_file2 character*128 output_file3 ! SJ33, May 2012 character*128 output_file4 ! SJ33, May 2012 character*16 out character*2 off,weak character*30 hkl character*12 hkl_int character*18 hkl_float character*24 datime integer ldt integer ntyp,i,j,numor,ii,iout,np radian=3.14159265/180.00 degrees=1./radian write(*,1000) 1000 format (/,' FIT CCL Files: Version 2.08, update Jun. 21, 2010', * //,' Data processing of files *.ccl from TriCS single' * ,' detector mode') call sys_setenv('DAT_DEFSPEC', 'TRICS') call dat_ask_filelist(input_files, ' ') output_file1=input_files i=index(output_file1, '.') if (i .eq. 0) i=index(output_file1, '[') if (i .eq. 0) i=index(output_file1, ',') if (i .eq. 0) i=index(output_file1, ' ') do j=i-1,1,-1 if (output_file1(j:j) .eq. '/') then goto 1002 endif enddo j=0 1002 if (i .gt. j+1) then output_file2=output_file1(j+1:i-1)//'_rafin.dat' output_file1=output_file1(j+1:i-1)//'.col' output_file3=output_file1(j+1:i-1)//'.comm' output_file4=output_file1(j+1:i-1)//'.incomm' else output_file2='rafin.dat' output_file1='upals.col' endif c write(*,1004) 1004 format (/,' function type to be used (0=gauss[def],6=strange)', * ' typ = ',$) read (5,'(i1)') ntyp if (ntyp .eq. 6) then write(*,'(/,'' width factor [1.0] = '',$)') read(5,'(f20.0)') wmult if (wmult .eq. 0.0) wmult=1.0 else wmult=1.0 ntyp=0 endif write(*,1013) 1013 format (/,' Monitor to be standardized , monitor=',$) read (5,'(f20.0)') yoverall if (yoverall.le.1.1) yoverall = 100000. write(*,1015) 1015 format (/, * ' plots (answer=n [no=def], y [yes,terminal], h [hard] ',$) read (5,'(a1)') answer plot = (answer.eq.'y' .or. answer .eq.'h') if (plot) then hard='y' if (answer.eq.'y') hard='n' endif write(*,1019) 1019 format (' Flat Background (y/n, default=y) : ',$) read (5,'(a1)')answer flat_backgr = (answer.ne.'n' .and. answer .ne.'N') if (ntyp .eq. 0) then write(*,1017) 1017 format * (//,' output intensities: 0=experimental (def), 1=gaussian ',$) read (5,'(i1)') Iout iexp= (iout.eq.0) else iexp= .true. endif if (iexp) then out='experimental ' else out='gaussian ' end if write(*,1018) 1018 format (' Lorentz-Correction to be applied (0 = 4-circle-geometry(default), 1 = tilt-geometry, 2 = not applied) ?: ',$) read (5,'(I1)') lorentz ! SJ33, June 2010, tilt option added (3 options now, def=0) if (lorentz.ne.1 .and. lorentz.ne.2) lorentz=0 ! SJ33, June 2010 write(*,1011) input_files,ntyp,wmult, * output_file1,output_file2,out,yoverall,flat_backgr,lorentz 1011 format (/,' ',78('-'),/,' input from file(s) :',A30,/, * ' function type ',i3,' width factor ',f6.3,//, * ' output: for UPALS :',A50,/, * ' for RAFIN :',A50,/, * 1X,A16,' intensities used'/, * ' monitor: ',f12.1,/, * ' Flat Background (yes/no): ',L33,/, * ' Lorentz-Correction (options: 0 = 4-circle-mode,', * ' 1 = tilt-mode, 2 = none)',/, * ' applied option No. ',i1,/, * ' No. 0 Icor=Iobs*sin(2theta)',/, * ' No. 1 Icor=Iobs*sin(Gama)*cos(nu)',/, * ' ',78('-'),/,' hit return to start ',$) c read (5,'(i1)') ii C if (ii.eq.1) goto 504 c if (ntyp.ne.0 .and. ntyp.ne.6) ntyp = 0 c c this ends the input for the program, the rest goes automatic c ------------------------------------------------------------ if (plot) then call sys_setenv('CHOOSER_PAN','9') ! 9 plots per page end if c c loop over all data following: iall =0 iincom =0 icom =0 icenter =0 open (unit=1,file=output_file1,status='unknown') open (unit=2,file=output_file2,status='unknown') open (unit=3,file=output_file3,status='unknown') open (unit=4,file=output_file4,status='unknown') c write headers for fullprof c write (3,1027)9.99,0,1,1,0.00,0.00,0.00 c write (4,1028)9.99,0,1,1,0.00,0.00,0.00 c1027 format (' title',/,22H(I6,3I10,2F10.2,4F8.2), /,F7.4,i5,i6,/,i1,3f5.2) c1028 format (' title',/,24H(I6,3F10.4,2F10.2,4F8.2),/,F7.4,I5,i6,/,i1,3f5.2) next_flag=0 ! start flag i=0 call fit_init c c ------------------------------------------------------ start loop -- 1500 continue call fit_dat_next(input_files, next_flag) if (next_flag .eq.0) goto 2222 ! last scan of last file read i=i+1 call fit_mon(yoverall) C set limits if peak to narrow: call fit_get_array('P', output, 20, ii) call fit_get_array('X', xx, np_max, np) c get information out of title: call fit_get_real('two_theta',th2) call fit_get_real('omega',omega_cal) call fit_get_real('chi',chi) call fit_get_real('phi',phi) call fit_get_real('h',h) call fit_get_real('k',k) call fit_get_real('l',l) call fit_get_real('temp',temperature) ! SJ33, Aug.30,1999 call fit_get_str('date', ldt, datime) cjs type *,' ***** temp=',temperature numor=i if (l.eq.88) then write(*,1007) 1007 format (' data numor not found, exit') goto 3000 endif c if (ntyp .eq. 6) then call fit_fun(6,1,wmult,0.0) ! select strange else call fit_fun(ntyp,0,0.0,0.0) ! select other (guassian) endif if (flat_backgr) call fit_set(2,0.,0.,0.,0.) ! flat background c range = scan with: range=xx(np)-xx(1) cjs type *,'******',output(5),output(6) if (ntyp .ne. 6) then if (output(6) .lt. 5.0*range/np) then c set position to center: call fit_set(3, (xx(np)+xx(1))/2.0, range/np, & xx(1)+range/6,xx(np)-range/6) c set fwhm fwhm = max(0.2,range/5) c arbitrary limits: minimal 3 steps, maximal 1 /3 range call fit_set(6,fwhm,fwhm/10.,3*range/np, range/3) else fwhm=max(output(6),5*range/np) call fit_set(6,fwhm,range/np 1 ,3*range/np,range/3) endif endif call fit_fit(0) c if (plot) call fit_plot(hard) ! hardcopy / terminal c call fit_get_array('P', output, 20, ii) call fit_get_array('E', output2, 20, ii) cjs type *,' output ',output cjs type *,' output2',output2 c c c get results now: if (ntyp.eq.0) then c function 1 (gauss: use only background, intensity=sum) omega_fit = output (3) omega_err = output2(3) off = ' ' if (iexp) then yint = output(8) dint = output2(8) else c modified 14.7.99 sj33 yint = output(5) dint = output2(5) endif bgr=output(1) fwhm=output(6) else c function 6 (strange) yint = output(7) dint = output2(7) omega_fit = output(3) omega_err = output2(3) bgr=output(1) fwhm=output(4) end if c c make notes for easy data inspections: delta = abs (omega_cal - omega_fit) if (delta.ge.0.3) off =' *' if (delta.ge.0.5) off ='**' weak = ' ' c c write on output files now: lor_cor = 1.00000 c L = 1 / SIN (2-Theta), page 156, Schwarzenbach, EPFL: c L applies to the calculated value of the intensities c Ical*L = Iobs c (measured intensity is always bigger than corrected one) if (lorentz.eq.0) lor_cor = sin (abs(th2)*radian) c c June 2010: SJ33 c tilt angle nu is read-in as chi from file *.col c correction formula: not used (3.59) on page 154, Schwarzenbach/Chapuis Cristallographie, Presses EPFL ! SJ33, June 2010 c used: Arndt and Wills, L=sin(gamma)*cos(nu), where gamma projection 2-theta, nu tilt yy1=lor_cor if (lorentz.eq.1) lor_cor = * sin(abs(th2)*radian)*cos(abs(chi)*radian) c * sqrt( (sin(abs(th2)*radian))**2 -(sin(abs(chi)*radian))**2 ) ! SJ33, June 2010 c write (6,*)' chi=',chi,th2,radian,sin(chi*radian),sin(th2*radian) c write (6,*)' Lorentz=',yy1,lor_cor c c if (abs(nint(h)-h) .gt. abs(h)*1e-4 .or. & abs(nint(k)-k) .gt. abs(k)*1e-4 .or. & abs(nint(l)-l) .gt. abs(l)*1e-4) then write(hkl, '(3(x,f9.4))') h, k, l write(hkl_float,'(3f6.2)') h, k, l iincom=iincom+1 iall=iall+1 comens=1 ! SJ33-May 2012 else write(hkl, '(3(x,i9))') nint(h), nint(k), nint(l) write(hkl_int,'(3i4)') nint(h), nint(k), nint(l) iall=iall+1 icom=icom+1 ! SJ33-May 2012 comens=0 endif c write (1, 1005) numor,hkl,lor_cor*yint,lor_cor*dint, & Th2/2.,omega_fit,chi,phi,temperature,yoverall/1000.,off,weak & ,bgr,fwhm,lor_cor,omega_err,datime(1:ldt) c if (comens.eq.0) then write (3, 1025) numor,hkl_int, lor_cor*yint,lor_cor*dint, * Th2/2,omega_fit,chi,phi ! SJ33-May 2012 else write (4, 1026) numor,hkl_float,lor_cor*yint,lor_cor*dint, * Th2/2,omega_fit,chi,phi endif ! SJ33-May 2012 1025 format (I6,A12,2F10.2,4f8.2) ! SJ33, nuc 1026 format (I4,A18,2F10.2,4F8.2) ! SJ33 mag 1005 format (I6,a,1x,f9.2,1x,f9.3,f7.2,1x,f8.3,1x,f7.2,1x,f7.2,1x, & F8.3,1x,F6.0,1x,2A2,1x,f6.1,1x,f5.3,1x,f6.4,1x,f6.3,1x,a) if (dint.ge.1) then s = yint/dint if (s.le.3) weak='w ' ! weak if (s.le.1) weak='vw' ! very weak if (s.ge.15 .and. yint.ge.100) then ! write strong hkl to rafin.dat write (2,1009) hkl,Th2,omega_fit,chi,phi,.2,temperature icenter=icenter+1 end if c 1009 format (a,f8.2,f8.3,2f8.2,1x,f5.3,f9.3) end if C goto 1500 c -------------------------------------------------------- end loop -- C 2222 write(*,1006) 1006 format (' Bye bye, normal end of program.') C C ------------------------------------------- C error exits here, does not pass label 2222: 3000 continue close (1) close (2) close (3) ! SJ33-May 2012 close (4) ! SJ33-May 2012 c c writing the end message for the user write (6,2040) write (6,2042) input_files(1:60) write (6,2041) iall,output_file1(1:30),icenter,output_file2(1:30), * icom,output_file3(1:30),iincom,output_file4 2040 format(///,' trics_ccl prepared the following files for YOU:',//) 2042 format(' TriCS-input-file used:',A60,//) 2041 format(i6,' reflections ',a30,' for JANA2006',/ * ,i6,' positions ',a30,' for rafin (to improve your UB using the dataset)',/ * ,I6,' reflections ',a30,' for Fullprof - commensurate data',/ * ,I6,' reflections ',a30,' for Fullprof - incommensurate data',//) c call fit_exit c end