Subroutine WRPOWDHIST for program GSAS2CIF
This subroutine is used to write parameters and results for each
powder histogram to the output CIF file.
See the gsas2cif documentation for an
explanation of this code.
SUBROUTINE WRPOWDHIST(IUCIF,IUEXP,IUTRM,IHST,HTYP,IUPRF,
1 LAM2,DAYTIME,ONEBLOCK,EXPRNAME,SAUTHOR)
C write the obs & calc powder diffractogram & reflection list
C************************************************************************
C
!PSEUDOCODE:
!CALLING ARGUMENTS:
INTEGER*4 IHST,IUCIF,IUEXP,IUTRM
INTEGER*4 IUPRF
CHARACTER*4 HTYP
REAL*4 LAM2 !Alpha_2 lambda
CHARACTER*20 DAYTIME
LOGICAL*4 ONEBLOCK ! true if the CIF will have one block
CHARACTER*20 EXPRNAME !Experiment name = name of data block
CHARACTER*24 SAUTHOR !A shortened version of the author name
!INCLUDE STATEMENTS:
INCLUDE '../INCLDS/ARRAYSZE.FOR'
INCLUDE '../INCLDS/SPGCOMI.FOR'
INCLUDE '../INCLDS/DISAGLCM.FOR'
!LOCAL VARIABLES:
INTEGER*4 NPHAS(9) !Phase existance flags
INTEGER*4 NPHASES ! number of phases in histogram
REAL*4 LAM1,ZERO,POLA,DIFC,DIFA
LOGICAL*4 MOREOBS ! true if there are more OBS points than calc
LOGICAL*4 FIXEDSTEP ! true for fixed step data
LOGICAL*4 FIXEDBKG ! true if fixed background points are used
LOGICAL*4 NEEDESD ! true if the ESD's are not SQRT(I)
REAL*4 DYDBK(99)
CHARACTER*80 TEXT
CHARACTER*80 TEXT1
CHARACTER*8 TYP
CHARACTER*20 BUFFER(10)
INTEGER*4 BUFLEN(10)
REAL*4 VALUE !General use value
REAL*4 FIRSTPT,LASTPT !Data range
REAL*4 STEPMIN,STEPMAX,STEP !Min, Max & avg step size
INTEGER*4 OFFSET !No. channels to be omitted at start of profile
INTEGER*4 ICLMP !Data compression factor
INTEGER*4 CHEKHST !Check sum of this histogram
INTEGER*4 ISAMP !Data samping factor
INTEGER*4 BAKTYP !Background function number
INTEGER*4 NUMBAK !Background number of terms
REAL*4 CONC(MAXELEM)
INTEGER*4 JMLT(MAXATM) !Atom site multiplicities
CHARACTER*12 KEYVAL !ISAM key
REAL*4 ATWT !Atomic weight
REAL*4 BLEN !Neutron scattering length
REAL*4 FFAC(9) !Xray form factor
REAL*4 FFAN(2,5) !Xray dispersion terms
REAL*4 ABSCO(7) !Absorption coefficients
REAL*4 MFAC(9) ! Neutron magnetic form factor
REAL*4 NFAC(9) ! Neutron magnetic form factor
REAL*4 BACKCOF(MAXBAK) !background coeffients
INTEGER*4 PRETRM !No. terms before diffuse terms in #9
INTEGER*4 TRMTYP(12) !Diffuse term types
CHARACTER*1 IAB !Absorption refinement flag
REAL*4 PHKL(3),RATIO,FRAC !M-D pref. orient.
REAL*4 COFF(MAXODF) !Spherical harmonic coefficients
INTEGER*4 INDX(3,MAXODF) !Sph. harmonic index
INTEGER*4 ISAMSYM !Sample symmetry number (1-4)
REAL*4 PAXIS(3) !Aniso. broadening axis
CHARACTER*1 SAXIS !=Y if stacking fault model is needed
REAL*4 UAXIS(3),VAXIS(3) !Stacking fault subcell vectors
INTEGER*4 PTYP,NPRF,PTA,PTB !Profile type & no. of coefficients
REAL*4 PCOF(36) !Profile coefficients
REAL*4 CTOF !Peak cutoff
!FUNCTION DEFINITIONS:
INTEGER*4 READEXP !ISAM file read function
CHARACTER*6 HSTKEY !ISAM key building routine
CHARACTER*6 HAPKEY !ISAM key building routine
CHARACTER*6 CRSKEY !ISAM key building routine
INTEGER*4 READPRF
LOGICAL*4 BTEST
INTEGER*4 LENCH !LENGTH of a character string
call OPNPRF(IUEXP,IHST,NCHAN,'SHARED',.FALSE.,IUPRF)
ISAM = READEXP(IUEXP,HSTKEY(IHST)//' NPHAS',TEXT)
READ (TEXT,'(9I5)') NPHAS
NPHASES = 0
do i=1,9
IF ( NPHAS(I).NE.0 ) THEN
NPHASES = NPHASES + 1
IPHAS = I
END IF
END DO
C======================================================================
C prepare for March-Dollase Preferred Orientation correction
C do any phases have a M-D correction?
IMD = 0
DO I = 1,9
IF ( NPHAS(I).NE.0 ) THEN
KEYVAL = HAPKEY(I,IHST)//'NAXIS '
IERR = READEXP(IUEXP,KEYVAL,TEXT)
READ(TEXT,'(I5)') NAXIS
JAX = 0
NUMF = 0
DO IAX=1,NAXIS
JAX = JAX+1
KEYVAL = HAPKEY(I,IHST)//'PREFO'//CHAR(48+JAX)
IERR = READEXP(IUEXP,KEYVAL,TEXT)
IF ( IERR.EQ.0 ) THEN
READ (TEXT,'(5F10.0)') RATIO,FRAC,(PHKL(K),K=1,3)
IF (RATIO .NE. 1) IMD = 1
END IF
END DO
END IF
END DO
C in the single-block case, need to also check for a spherical harmonic
IODF = 0
IF (ONEBLOCK) THEN
DO I = 1,9
IF ( NPHAS(I).NE.0 ) THEN
KEYVAL = CRSKEY(I)//'ODF '
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
READ (TEXT,'(2I5)') NORD,NODFCOF
IF (NODFCOF .GT. 0) IODF = 1
END IF
END DO
END IF
C======================================================================
IF (.NOT. ONEBLOCK) THEN
WRITE(IUCIF,'(A)') '# phase table'
WRITE(IUCIF,'(A)') 'loop_ _pd_phase_id'
WRITE(IUCIF,'(10X,A)') '_pd_phase_block_id'
WRITE(IUCIF,'(10x,A)') '_pd_phase_mass_%'
IF (IMD .GT. 0) WRITE(IUCIF,'(10X,A)')
1 '_pd_proc_ls_pref_orient_corr'
WRITE(IUCIF,'(10X,A)') '_pd_proc_ls_profile_function'
WRITE(IUCIF,'(10X,A)') '_pd_proc_ls_peak_cutoff'
DO I=1,9
IF ( NPHAS(I).NE.0 ) THEN
C_pd_phase_block_id
WRITE(IUCIF,'(2X,I1,2X,2A,I1,4A)') I,
1 DAYTIME(1:16)//'|',
1 EXPRNAME(1:LENCH(EXPRNAME))//'_phase',I,'|',
1 SAUTHOR(:LENCH(SAUTHOR)),'||'
C_pd_phase_mass_% (from phase fractions)
TEXT = ' '
KEYVAL = HAPKEY(I,IHST)//'MASSFR'
IERR = READEXP(IUEXP,KEYVAL,TEXT)
IF (TEXT .NE. ' ') THEN
READ(TEXT,'(2F10.4)') WTFR,SIGW
CALL FESD(WTFR*100., SIGW*100., TEXT, LN)
WRITE (IUCIF,'(10x,A)') TEXT(:LN)
ELSE
WRITE (IUCIF,'(10x,A)') '?'
ENDIF
C_pd_proc_ls_pref_orient_corr
IF (IMD .GT. 0) THEN
WRITE(IUCIF,'(2A)') ';',' March-Dollase'
TEXT = ' '
KEYVAL = HAPKEY(I,IHST)//'NAXIS '
IERR = READEXP(IUEXP,KEYVAL,TEXT)
READ(TEXT,'(I5)') NAXIS
JAX = 0
NUMF = 0
DO IAX=1,NAXIS
JAX = JAX+1
KEYVAL = HAPKEY(I,IHST)//'PREFO'//CHAR(48+JAX)
IERR = READEXP(IUEXP,KEYVAL,TEXT)
IF ( IERR.EQ.0 ) THEN
READ (TEXT,'(5F10.0)') RATIO,FRAC,
1 (PHKL(K),K=1,3)
IF (NAXIS .EQ. 1) THEN
WRITE(IUCIF,'(A,I2,A,F10.5,3(2x,A,F6.3))')
1 ' AXIS ',IAX,' Ratio=',RATIO,
1 'h=',PHKL(1),'k=',PHKL(2),'l=',PHKL(3)
ELSE
WRITE(IUCIF,'(A,I2,2(A,F10.3),3(2x,A,F6.3))')
1 ' AXIS ',IAX,
1 ' Ratio=',RATIO,' Frac',FRAC,
1 'h=',PHKL(1),'k=',PHKL(2),'l=',PHKL(3)
END IF
END IF
END DO
WRITE(IUCIF,'(A))') ';'
END IF
C_pd_proc_ls_profile_function
WRITE(IUCIF,'(A))') ';'
CALL SYMMINP(IUEXP,I,ICEN,NSYM,LCENT,NCV,LAUE,NPOL,
1 NAXIS,JRT,CEN,SPG)
KEYVAL = HAPKEY(I,IHST)//'PRCF '
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
READ(TEXT,'(2I5,F10.5,I5)') PTYP,NPRF,CTOF,IDAMP
DO K=1,36
PCOF(K) = 0.0
END DO
NREC = (MPRF-1)/4+1
DO IREC=1,NREC
WRITE(KEYVAL(12:12),'(I1)') IREC
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
IBEG = (IREC-1)*4+1
READ(TEXT,'(4E15.6)') (PCOF(K),K=IBEG,IBEG+3)
END DO
C now for some pain... list the profile terms for all of Bob's masterpieces
CALL LISTPRF(IUCIF,NPRF,PTYP,PCOF,LAUE,NAXIS,HTYP,CTOF)
KEYVAL = CRSKEY(I)//'SPAXIS'
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
IF ( ISAM.EQ.0 ) THEN
READ(TEXT,'(3F5.0,4X,A,6F5.0)') PAXIS,SAXIS,UAXIS,VAXIS
IF ( SAXIS.EQ.'Y' ) THEN
WRITE(IUCIF,3) PAXIS,UAXIS,VAXIS
3 FORMAT(' Stacking fault sublattice vectors:',/,
1 10x,2(3F6.1,';'),3F6.1)
ELSE
WRITE(IUCIF,'(A,3F6.1)')
1 ' Aniso. broadening axis',PAXIS
END IF
END IF
WRITE(IUCIF,'(A)') ';'
C_pd_proc_ls_peak_cutoff
WRITE(IUCIF,'(10X,F8.5)') CTOF
END IF
END DO
END IF
C======================================================================
C loop over atom types and report the scattering factor info
C include atom amounts, if one histogram and one phase (note that phase info
C was loaded in WRITEPHASE)
IF (ONEBLOCK) THEN
C determine unit cell contents
ISAM = READEXP(IUEXP,' EXPR NATYP',TEXT)
READ (TEXT,'(9I5)') NELEM
DO I=1,NELEM
conc(i) = 0.0
END DO
DO I=1,NATOM
CALL SYTSYM (XYZ(1,I),ICEN,NSYM,JRT,NCV,CEN,LAUE, !Get site symmetries and multiplicities
1 JMLT(I),JSYM,IOPRTNS,STSYM(I))
J = ID(I) !Get the atom type count flag
conc(J) = conc(j) + JMLT(I)*FRACT(I)
END DO
END IF
C loop header
WRITE(IUCIF,'(/a5,1x,A)') 'loop_', '_atom_type_symbol'
IF (ONEBLOCK) THEN
WRITE(IUCIF,'(6x,A)') '_atom_type_number_in_cell'
END IF
IRAD = -1
IF (HTYP(2:2) .eq. 'X') THEN
WRITE(IUCIF,'(6x,A)') '_atom_type_scat_dispersion_real'
WRITE(IUCIF,'(6x,A)') '_atom_type_scat_dispersion_imag'
WRITE(IUCIF,'(6x,''_atom_type_scat_Cromer_Mann_'',A)') 'a1',
1 'a2','a3', 'a4', 'b1', 'b2', 'b3', 'b4', 'c'
ISAM = READEXP(IUEXP,HSTKEY(IHST)//' IRAD ',TEXT)
READ (TEXT,'(I5)') IRAD
ELSEIF (HTYP(2:2) .eq. 'N') THEN
WRITE(IUCIF,'(6x,A)') '_atom_type_scat_length_neutron'
END IF
WRITE(IUCIF,'(6x,A)') '_atom_type_scat_source'
C _atom_type_description
ISAM = READEXP(IUEXP,' EXPR NATYP',TEXT)
READ (TEXT,'(9I5)') NELEM
IF (NELEM .LE. 0) THEN
PRINT '(A)','Error -- This experiment contains no atom types.'
STOP 'EXPR NATYP Error'
END IF
DO j=1,NELEM
CALL VSTRNG(ATYPE(j),LENCH(ATYPE(j)),.true.,.false.)
text = ATYPE(j)(1:8)
ln = 10
IF (ONEBLOCK) THEN
CALL FESD(conc(j), -0.01, text(ln:), ln)
ln = ln + 2
END IF
CALL RDTYPDT(IUEXP,ATYPE(j),ATWT,BLEN,FFAC,FFAN,ABSCO,MFAC,
1 NFAC,GFAC)
IF (HTYP(2:2) .eq. 'X') THEN
IF (IRAD .GE. 1 .and. IRAD .LE. 5) THEN
write(IUCIF,'(2x,a,2f9.3)') text(:LENCH(text)),
1 (FFAN(I,IRAD),I=1,2)
ELSE IF (IRAD .EQ. 0) THEN
FF1 = 0
FF2 = 0
IREC = 0
ISAM = 0
C loop through the anomolous f' & f'' for values
DO WHILE (ISAM .EQ. 0 .AND. IREC .LT. 9)
IREC = IREC + 1
KEYVAL = HSTKEY(IHST)//'FFANS '
WRITE(KEYVAL(12:12),'(I1)') IREC
ISAM = READEXP(IUEXP,KEYVAL,TEXT1)
IF ( ISAM.EQ.0 ) THEN
READ(TEXT1,'(2X,A8,2F10.0)') TYP,FF1A,FF2A
IF (TYP .EQ. ATYPE(j)) THEN
FF1 = FF1A
FF2 = FF2A
ISAM = -1
END IF
END IF
END DO
write(IUCIF,'(2x,a,2f9.3)') text(:LENCH(text)),FF1,FF2
END IF
C now write the coeff. for the scattering curve
LN = 1
TEXT = ' '
DO I=1,9
IF (FFAC(I) .GE. 10000 .OR. FFAC(I) .LE. -1000) THEN
WRITE(TEXT(LN:),'(F8.1)') FFAC(I)
ELSEIF (FFAC(I) .GE. 1000 .OR. FFAC(I) .LE. -100) THEN
WRITE(TEXT(LN:),'(F8.2)') FFAC(I)
ELSEIF (FFAC(I) .GE. 100 .OR. FFAC(I) .LE. -10) THEN
WRITE(TEXT(LN:),'(F8.3)') FFAC(I)
ELSEIF (FFAC(I) .GE. 10 .OR. FFAC(I) .LE. 0) THEN
WRITE(TEXT(LN:),'(F8.4)') FFAC(I)
ELSE
WRITE(TEXT(LN:),'(F8.5)') FFAC(I)
ENDIF
LN = LN + 8
END DO
WRITE(IUCIF,'(A)') TEXT(:LN)
ELSE IF (HTYP(2:2) .eq. 'N') THEN
write(IUCIF,'(2x,a,1x,F8.4)') text(:LENCH(text)),BLEN
END IF
C_atom_type_scat_source
write(IUCIF,'(2x,a)') 'International_Tables_Vol_C'
END DO
C======================================================================
IF (HTYP(2:2) .eq. 'X') THEN
CALL WRVAL(IUCIF,'_diffrn_radiation_probe','x-ray')
ELSE IF (HTYP(2:2) .eq. 'N') THEN
CALL WRVAL(IUCIF,'_diffrn_radiation_probe','neutron')
ELSE
PRINT '(A)','Unexpected data type for _diffrn_radiation'//
1 '_probe histogram #',IHST
CALL WRVAL(IUCIF,'_diffrn_radiation_probe','?')
END IF
ISAM = readexp(IUEXP, HSTKEY(IHST)//' CHANS',text)
READ(TEXT,'(20x,I10,10x,i10)') nchans,mchans
ISAM = READEXP(IUEXP, HSTKEY(IHST)//' ICONS',TEXT)
IF (HTYP(3:3) .eq. 'T') THEN
READ(TEXT,'(3F10.0)') DIFC,DIFA,ZERO
ELSE IF (HTYP(2:2) .eq. 'N') THEN
READ (TEXT,'(3f10.0,25x,f10.0)') lam1,lam2,zero,ratio
CALL FESD(LAM1, -.00001, text, ln)
CALL WRVAL(IUCIF, '_diffrn_radiation_wavelength' ,text)
ELSE
READ (TEXT,'(3f10.0,10x,f10.0,i5,f10.0)') lam1,lam2,zero,
1 pola,ipola,ratio
C Bob -- I don't know how to treat case of IPOLA != 0
IF (IPOLA .EQ. 0) THEN
CALL FESD(POLA, -.01, text, ln)
CALL WRVAL(IUCIF, '_diffrn_radiation_polarisn_ratio' ,text)
ELSE
CALL WRVAL(IUCIF, '_diffrn_radiation_polarisn_ratio' ,'?')
END IF
IF (LAM2 .eq. 0.0) then
CALL FESD(LAM1, -.00001, text, ln)
CALL WRVAL(IUCIF, '_diffrn_radiation_wavelength' ,text)
CALL WRVAL(IUCIF, '_diffrn_radiation_type', '?')
ELSE
C always assume Ka1 & Ka2 if two wavelengths are present
WRITE(IUCIF,'(/a)') 'loop_'
WRITE(IUCIF,'(6x,A)') '_diffrn_radiation_wavelength'
WRITE(IUCIF,'(6x,A)') '_diffrn_radiation_wavelength_wt'
WRITE(IUCIF,'(6x,A)') '_diffrn_radiation_type'
WRITE(IUCIF,'(6x,A)') '_diffrn_radiation_wavelength_id'
WRITE(IUCIF,'(20x,f10.6,5x,f6.3,3x,A,i3)')
1 LAM1,1.0,'K\\a~1~',1,
1 LAM2,ratio,'K\\a~2~',2
END IF
END IF
TEXT = ' ? ?'
ISAM = READEXP(IUEXP,HSTKEY(IHST)//' RPOWD',TEXT)
CALL WRVAL(IUCIF,'_pd_proc_ls_prof_R_factor',TEXT(11:20))
CALL WRVAL(IUCIF,'_pd_proc_ls_prof_wR_factor',TEXT(1:10))
TEXT = ' ?'
ISAM = READEXP(IUEXP,HSTKEY(IHST)//' WREXP',TEXT)
CALL WRVAL(IUCIF,'_pd_proc_ls_prof_wR_expected',TEXT(1:10))
TEXT = ' ?'
ISAM = READEXP(IUEXP,HSTKEY(IHST)//' R-FAC',TEXT)
CALL WRVAL(IUCIF,'_refine_ls_R_Fsqd_factor',TEXT(6:15))
C document the background function used
KEYVAL = HSTKEY(IHST)//' TRNGE'
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
IF ( ISAM.NE.0 ) TEXT = ' 1.0 100.0'
READ(TEXT,'(2F10.0)') TTMIN,TTMAX
KEYVAL = HSTKEY(IHST)//'BAKGD '
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
IF (ISAM .EQ. 0) THEN
READ(TEXT,'(2I5,15X,I5,12I1)') BAKTYP,NUMBAK,PRETRM,TRMTYP
WRITE(IUCIF,'(/a)') '_pd_proc_ls_background_function'
WRITE(IUCIF,'(2a,I2,a,I3,a)') ';',
1 ' GSAS Background function number',BAKTYP,
1 ' with',NUMBAK,' terms.'
IF ( BAKTYP.EQ.1 ) THEN
WRITE(IUCIF,'(A)') ' Shifted Chebyshev function of 1st kind'
ELSE IF ( BAKTYP.EQ.2 ) THEN
WRITE(IUCIF,'(A)') ' Cosine Fourier series'
ELSE IF ( BAKTYP.EQ.3 ) THEN
WRITE(IUCIF,'(A)') ' Real space distribution function'
ELSE IF ( BAKTYP.EQ.4 ) THEN
WRITE(IUCIF,'(A)') ' Power series in Q**2n/n!'
ELSE IF ( BAKTYP.EQ.5 ) THEN
WRITE(IUCIF,'(A)') ' Power series in n!/Q**2n'
ELSE IF ( BAKTYP.EQ.6 ) THEN
WRITE(IUCIF,'(A)') ' Power series in Q**2n/n! and n!/Q**2n'
ELSE IF ( BAKTYP.EQ.7 ) THEN
WRITE(IUCIF,'(A)') ' Linear interpolation'
ELSE IF ( BAKTYP.EQ.8 ) THEN
WRITE(IUCIF,'(A)') ' Reciprocal interpolation'
ELSE IF ( BAKTYP.EQ.9 ) THEN
WRITE(IUCIF,'(A)') ' Diffuse scattering function'
END IF
NREC = (NUMBAK-1)/4+1
IBEG = 1
DO IREC=1,NREC
WRITE(KEYVAL(12:12),'(I1)')IREC
IFIN = MIN(IBEG+3,NUMBAK)
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
READ (TEXT,'(4E15.6)') (BACKCOF(I),I=IBEG,IFIN)
WRITE(IUCIF,'(5x,4(I2,A,1PG15.6))')
1 (I,':',BACKCOF(I),I=IBEG,IFIN)
IBEG = IBEG+4
END DO
WRITE(IUCIF,'(a)') ';'
END IF
C handle absorption/roughness correction
KEYVAL = HSTKEY(IHST)//'ABSCOR'
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
IF (ISAM .EQ. 0 .AND. TEXT(20:20).EQ.'.' ) THEN !Ignore old record format
READ(TEXT,'(2E15.6,4X,A,2I5)') ABSCOR1,ABSCOR2,IAB,
1 IDAMP,IABTYP
WRITE(IUCIF,'(/a)') '_exptl_absorpt_process_details'
WRITE(IUCIF,'(2a,I2)') '; GSAS Absorption/surface roughness',
1 ' correction: function number',IABTYP
IF ( IABTYP.EQ.0 .AND. ABSCOR1 .EQ. 0) THEN
WRITE(IUCIF,'(A)') ' No correction is applied.'
ELSE IF ( IABTYP.EQ.0 ) THEN
WRITE(IUCIF,'(A)') ' Debye-Scherrer absorption correction'
WRITE(IUCIF,'(A,G15.5)') 'Term (= MU.r/wave) = ',ABSCOR1
ELSE IF ( IABTYP.EQ.1 ) THEN
WRITE(IUCIF,'(A)') ' Linear absorption correction'
WRITE(IUCIF,'(A,2G15.5)') 'Term = ',ABSCOR1
ELSE IF ( IABTYP.EQ.2 ) THEN
WRITE(IUCIF,'(A)') ' Surface roughness abs. correction'//
1 ' (Pitschke, et al.)'
WRITE(IUCIF,'(A,2G15.5)') 'Terms = ',ABSCOR1,ABSCOR2
ELSE IF ( IABTYP.EQ.3 ) THEN
WRITE(IUCIF,'(A)') ' Surface roughness abs. correction'//
1 ' (Suortti)'
WRITE(IUCIF,'(A,2G15.5)') 'Terms = ',ABSCOR1,ABSCOR2
ELSE IF ( IABTYP.EQ.4 ) THEN
WRITE(IUCIF,'(A)') ' Flat plate transmission absorption'//
1 ' correction'
WRITE(IUCIF,'(A,2G15.5)') 'Terms = ',ABSCOR1,ABSCOR2
END IF
IF (IAB .EQ. 'Y') THEN
WRITE(IUCIF,'(A,2G15.5)') 'Correction is refined.'
ELSE IF (IABTYP.NE.0 .OR. ABSCOR1 .NE. 0) THEN
WRITE(IUCIF,'(A,2G15.5)') 'Correction is not refined.'
END IF
WRITE(IUCIF,'(a)') ';'
END IF
C probably not needed
C _exptl_absorpt_correction_type 'shelx76 gaussian'
C not exactly appropriate
C _exptl_absorpt_coefficient_mu 0.59 (unknown in GSAS? -- empirical?)
C _pd_char_atten_coef_mu_obs
C show range of applied corrections
C absorption
TEXT = ' '
ISAM = READEXP(IUEXP, HSTKEY(IHST)//'TRMNMX',TEXT)
IF (TEXT .NE. ' ') THEN
CALL WRVAL(IUCIF,'_exptl_absorpt_correction_T_min',TEXT(1:10))
CALL WRVAL(IUCIF,'_exptl_absorpt_correction_T_max',TEXT(11:20))
ELSE
CALL WRVAL(IUCIF,'_exptl_absorpt_correction_T_min','?')
CALL WRVAL(IUCIF,'_exptl_absorpt_correction_T_max','?')
ENDIF
C extinction
TEXT = ' '
ISAM = READEXP(IUEXP, HSTKEY(IHST)//'EXMNMX',TEXT)
IF (TEXT .NE. ' ') THEN
WRITE(IUCIF,'(A)') '# Extinction correction'
CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_min',TEXT(1:10))
CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_max',TEXT(11:20))
ENDIF
IF (ONEBLOCK .AND. IMD+IODF .GT. 0) THEN
I = IPHAS
WRITE(IUCIF,'(A)') '_pd_proc_ls_pref_orient_corr'
IF (IMD .GT. 0) THEN
WRITE(IUCIF,'(2A)') ';',' March-Dollase'
ELSE
WRITE(IUCIF,'(2A)') ';',' Spherical Harmonic ODF'
END IF
IF (IMD .GT. 0) THEN
TEXT = ' '
KEYVAL = HAPKEY(I,IHST)//'NAXIS '
IERR = READEXP(IUEXP,KEYVAL,TEXT)
READ(TEXT,'(I5)') NAXIS
JAX = 0
NUMF = 0
DO IAX=1,NAXIS
JAX = JAX+1
KEYVAL = HAPKEY(I,IHST)//'PREFO'//CHAR(48+JAX)
IERR = READEXP(IUEXP,KEYVAL,TEXT)
IF ( IERR.EQ.0 ) THEN
READ (TEXT,'(5F10.0)') RATIO,FRAC,
1 (PHKL(K),K=1,3)
IF (NAXIS .EQ. 1) THEN
WRITE(IUCIF,'(A,I2,A,F10.5,3(2x,A,F6.3))')
1 ' AXIS ',IAX,' Ratio=',RATIO,
1 'h=',PHKL(1),'k=',PHKL(2),'l=',PHKL(3)
ELSE
WRITE(IUCIF,'(A,I2,2(A,F10.3),3(2x,A,F6.3))')
1 ' AXIS ',IAX,
1 ' Ratio=',RATIO,' Frac',FRAC,
1 'h=',PHKL(1),'k=',PHKL(2),'l=',PHKL(3)
END IF
END IF
END DO
END IF
IF (IODF .GT. 0 .and. IMD .NE. 0) THEN
WRITE(IUCIF,'(/A)') ' **** Spherical Harmonic ODF ****'
END IF
IF (IODF .GT. 0) THEN
KEYVAL = CRSKEY(I)//'ODF '
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
READ (TEXT,'(3I5)') NORD,NODFCOF,ISAMSYM
WRITE(IUCIF,'(A,I3)')
1 ' Spherical harmonic order=',NORD
IF ( ISAMSYM.EQ.1 ) THEN
WRITE(IUCIF,'(A)') ' No sample symmetry'
ELSE IF ( ISAMSYM.EQ.2 ) THEN
WRITE(IUCIF,'(2A)') ' The sample symmetry is:',
1 ' 2/m (shear texture)'
ELSE IF ( ISAMSYM.EQ.3 ) THEN
WRITE(IUCIF,'(2A)') ' The sample symmetry is:',
1 ' mmm (rolling texture)'
ELSE IF ( ISAMSYM.EQ.0 ) THEN
WRITE(IUCIF,'(2A)') ' The sample symmetry is:',
1 ' cylindrical (fiber texture)'
END IF
NREC = 0
IF ( NODFCOF.GT.0 ) NREC = (NODFCOF-1)/6+1
IBEG = 1
DO IREC=1,NREC
WRITE(KEYVAL(10:12),'(I2,A)')IREC,'A'
IFIN = MIN(IBEG+5,NODFCOF)
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
READ (TEXT,'(6(I4,2I3))') ((INDX(K,J),K=1,3),
1 J=IBEG,IFIN)
KEYVAL(12:12) = 'B'
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
READ (TEXT,'(6(F10.0))') (COFF(K),K=IBEG,IFIN)
IBEG = IBEG+6
END DO
DO J=1,NODFCOF
WRITE(IUCIF,'(A,3I3,3x,A,F10.4)')
1 ' Index =',(INDX(K,J),K=1,3),
1 'Coeff=',COFF(J)
END DO
END IF
TEXT = ' '
ISAM = READEXP(IUEXP, HSTKEY(IHST)//'ODMNMX',TEXT)
IF (TEXT .NE. ' ') THEN
WRITE(IUCIF,'(4A)')
1 ' Prefered orientation correction range: Min=',
1 TEXT(1:10),', Max=',TEXT(11:20)
ENDIF
WRITE(IUCIF,'(A)') ';'
END IF
C_pd_proc_ls_profile_function
IF (ONEBLOCK) THEN
I = IPHAS
WRITE(IUCIF,'(2(/,A))') '_pd_proc_ls_profile_function',';'
CALL SYMMINP(IUEXP,I,ICEN,NSYM,LCENT,NCV,LAUE,NPOL,
1 NAXIS,JRT,CEN,SPG)
KEYVAL = HAPKEY(I,IHST)//'PRCF '
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
READ(TEXT,'(2I5,F10.5,I5)') PTYP,NPRF,CTOF,IDAMP
DO K=1,36
PCOF(K) = 0.0
END DO
NREC = (MPRF-1)/4+1
DO IREC=1,NREC
WRITE(KEYVAL(12:12),'(I1)') IREC
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
IBEG = (IREC-1)*4+1
READ(TEXT,'(4E15.6)') (PCOF(K),K=IBEG,IBEG+3)
END DO
C now for some pain... list the profile terms for all of Bob's masterpieces
CALL LISTPRF(IUCIF,NPRF,PTYP,PCOF,LAUE,NAXIS,HTYP,CTOF)
KEYVAL = CRSKEY(I)//'SPAXIS'
ISAM = READEXP(IUEXP,KEYVAL,TEXT)
IF ( ISAM.EQ.0 ) THEN
READ(TEXT,'(3F5.0,4X,A,6F5.0)') PAXIS,SAXIS,UAXIS,VAXIS
IF ( SAXIS.EQ.'Y' ) THEN
WRITE(IUCIF,3) PAXIS,UAXIS,VAXIS
ELSE
WRITE(IUCIF,'(A,3F6.1)')
1 ' Aniso. broadening axis',PAXIS
END IF
END IF
WRITE(IUCIF,'(A)') ';'
WRITE(IUCIF,'(A,F8.5)') '_pd_proc_ls_peak_cutoff',CTOF
END IF
C use current time/date here
CALL WRVAL(IUCIF, '_pd_proc_info_datetime', daytime)
CALL WRVAL(IUCIF, '_pd_calc_method', 'Rietveld Refinement')
DO I=1,10
BUFLEN(I) = 0
END DO
C put the intensity data on a scratch file, so that we can find the length
C of each column; then we can line up numbers so they look pretty
CALL GETUNIT(IUSCRT)
OPEN(IUSCRT)
C is this time-of-flight?
IF (HTYP(3:3) .eq. 'T') THEN
FIXEDSTEP = .false. ! true for fixed step data
NEEDESD = .true.
ELSE
C check through the data to check if the step size is fixed and
C get the starting, ending angles & step angles/channel while doing this
J = 1
IREC = 0
STEPMIN = 0
STEPMAX = 0
LASTPT = -1
C are the intensity values scaled? Assume no & scan through the histogram
NEEDESD = .false.
do while (J .ne. 0)
IREC = IREC + 1
J = READPRF(IUPRF,IREC,ICODE,FIRSTPT,YO,YC,YI,YB,YW,CHWDT,
1 MIN1,MIN2)
END DO
K = 1
IF (YW .GT. 0 .AND.
1 (YO*YW .LT. .95 .OR. YO*YW .GT. 1.05))
1 NEEDESD = .true.
ISAM = readexp(IUEXP, HSTKEY(IHST)//' CHANS',text)
IF ( ISAM.EQ.0 ) READ(text,'(5I10,I5)') OFFSET,ICLMP,
1 NCHANS,CHEKHST,MCHANS,ISAMP
IF ( ISAMP.EQ.0 ) ISAMP = 1
DO I = IREC+1,NCHANS
J = READPRF(IUPRF,I,ICODE,T2,YO,YC,YI,YB,YW,CHWDT,MIN1,MIN2)
if (j .eq. 0) then
IF (YW .GT. 0 .AND.
1 (YO*YW .LT. .95 .OR. YO*YW .GT. 1.05))
1 NEEDESD = .true.
C Is this the second defined point?
IF (LASTPT .EQ. -1) THEN
STEPMIN = ABS(T2 - FIRSTPT)
STEPMAX = ABS(T2 - FIRSTPT)
ELSE
STEPMIN = MIN(STEPMIN,ABS(T2 - LASTPT))
STEPMAX = MIN(STEPMAX,ABS(T2 - LASTPT))
END IF
irec = I
LASTPT = t2
k = k + 1
END IF
END DO
C treat a <1% variation in stepsize as fixed step size
IF ( (STEPMAX-STEPMIN)/STEPMAX .GT. 0.01) THEN
FIXEDSTEP = .false.
ELSE
C round step to nearest .001 degree
STEP = NINT(100.*((LASTPT - FIRSTPT)/(k-1.)))/100.
FIXEDSTEP = .true.
END IF
END IF
C do we have any fixed background points?
FIXEDBKG = .false. ! true if fixed background points are used
ISAM = READEXP(IUEXP,HSTKEY(IHST)//'FXB 1',TEXT)
IF ( ISAM.NE.0 ) FIXEDBKG = .true.
C do we have a different number of experimental data points then in the
C refined histogram? Caused by data compression or sampling
C for now, ignore dropped points at the beginning of the diffraction pattern
IF (ICLMP .GT. 1 .OR. ISAMP .GT. 1) THEN
MOREOBS = .true. ! there are more OBS points than calc
ELSE
MOREOBS = .false. ! same number of OBS & CALC points
END IF
C**********************************************************************
C loop for raw data (if needed)
IF (MOREOBS) THEN
WRITE(IUCIF,'(/,A)') '#---- raw data loop -----'
CALL WRITERAWDATA(IUEXP,IUCIF,IHST,HTYP,FIXEDSTEP)
WRITE(IUCIF,'(/,A)') '#---- calculated data loop -----'
ELSE
WRITE(IUCIF,'(/,A)') '#---- raw/calc data loop -----'
END IF
C**********************************************************************
C loop over computed histogram & optionally list the observed as well
IF (FIXEDSTEP) THEN
C starting, ending angles & step for OBS & CALC -- N.B. Always 2theta
IF (.not. MOREOBS) then
CALL FESD(FIRSTPT/100., -abs(STEP/10000.), text,ln)
CALL WRVAL(IUCIF, '_pd_meas_2theta_range_min', text)
CALL FESD(LASTPT/100., -abs(STEP/10000.), text,ln)
CALL WRVAL(IUCIF, '_pd_meas_2theta_range_max', text)
CALL FESD(STEP/100., -abs(STEP/10000.), text,ln)
CALL WRVAL(IUCIF, '_pd_meas_2theta_range_inc', text)
END IF
C zero correct & convert from centidegrees
FIRSTPT = (FIRSTPT - zero)/100.
LASTPT = (LASTPT - zero)/100.
CALL FESD(FIRSTPT, -abs(STEP/10000.), text,ln)
CALL WRVAL(IUCIF, '_pd_proc_2theta_range_min', text)
CALL FESD(LASTPT, -abs(STEP/10000.), text,ln)
CALL WRVAL(IUCIF, '_pd_proc_2theta_range_max', text)
CALL FESD(STEP/100., -abs(STEP/10000.), text,ln)
CALL WRVAL(IUCIF, '_pd_proc_2theta_range_inc', text)
END IF
C**********************************************************************
C write the header for the main data loop
WRITE(IUCIF,'(/,A)') 'loop_'
IF (HTYP(3:3) .eq. 'T') THEN
IF (.not. MOREOBS) WRITE(IUCIF,'(6x,A)')
1 '_pd_meas_time_of_flight'
WRITE(IUCIF,'(6x,A)') '_pd_proc_d_spacing'
ELSE IF (.not. FIXEDSTEP) THEN
IF (.not. MOREOBS) WRITE(IUCIF,'(6x,A)') '_pd_meas_2theta_scan'
IF (MOREOBS .or. ZERO .ne. 0.)
1 WRITE(IUCIF,'(6x,A)') '_pd_proc_2theta_corrected'
END IF
C which intensity variable is needed?
IF (MOREOBS) THEN
WRITE(IUCIF,'(6x,A)') '_pd_proc_intensity_total'
NEEDESD = .TRUE. ! this requires SU's since _total is assumed to not be counts
ELSE IF (NEEDESD) THEN
WRITE(IUCIF,'(6x,A)') '_pd_meas_intensity_total'
ELSE
WRITE(IUCIF,'(6x,A)') '_pd_meas_counts_total'
END IF
WRITE(IUCIF,'(6x,A)') '_pd_proc_ls_weight'
WRITE(IUCIF,'(6x,A)') '_pd_proc_intensity_bkg_calc'
C for now ignore fixed background points
! IF (FIXBACK) WRITE(IUCIF,'(6x,A)') '_pd_proc_intensity_fix_bkg'
WRITE(IUCIF,'(6x,A)') '_pd_calc_intensity_total'
KEYVAL = HSTKEY(IHST)
IF ( HTYP(3:3).EQ.'T' ) THEN !TOF neutron data
ISAM = READEXP(IUEXP,KEYVAL(1:6)//'BNKPAR',TEXT)
READ(TEXT,'(10X,F10.0)') WAVE
DIFC1000 = DIFC/1000.
ISAM = READEXP(IUEXP,KEYVAL(1:6)//'I ITYP',TEXT)
READ(TEXT,'(15X,F10.4)') TMAX
TMAX = 180.0/TMAX
ELSE IF ( HTYP(3:3).EQ.'C' ) THEN
WAVE = LAM1
TMAX = 180.0
ELSE IF ( HTYP(3:3).EQ.'E' ) THEN
WAVE = LAM1
TMAX = 250.0
END IF
C now loop through the data
npoint = 0
DO I=1,nchans
j = 0
k = READPRF(IUPRF,I,ICODE,T1,YO,YC,YI,YB,YW,CHWDT,MIN1,MIN2)
if (k .eq. 0) then
C compute the background
npoint = npoint + 1
IF ( HTYP(3:3).EQ.'T' ) THEN ! TOF
T1A = T1/1000.
ELSE IF (HTYP(3:3).EQ.'C') THEN ! constant wavelength
T1A = T1/100.
ELSE IF ( HTYP(3:3).EQ.'E' ) THEN ! energy dispersive x-ray
T1A = T1
END IF
CALL CALCBAK(HTYP,TMAX,DIFC1000,WAVE,TTMIN,TTMAX,
1 BAKTYP,NUMBAK,BACKCOF,PRETRM,TRMTYP,T1A,YB1)
C add the fixed and computed background
YB1 = YB + YB1
C----------------------------------------------------------------------
C report the scan variable, if TOF or not fixed step
IF (HTYP(3:3) .eq. 'T') THEN
C _pd_meas_time_of_flight
IF (.not. MOREOBS) THEN
J = J + 1
CALL FESD(t1, -0.01, buffer(j),ln) !The LANSCE T-O-F clocks run @ 20M Hz
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
END IF
C _pd_proc_d_spacing
TMP = TOFTOD(HTYP,DIFC,DIFA,ZERO,1.0,T1,VALUE)
T1 = TMP
J = J + 1
LN = -1
CALL FESD(t1, -t1*.0001, buffer(j),ln)
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
ELSE IF (.not. FIXEDSTEP) THEN
C _pd_meas_2theta_scan
IF (.not. MOREOBS) THEN
J = J + 1
LN = -1
CALL FESD(t1/100., -0.001, buffer(j),ln)
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
END IF
C _pd_proc_2theta_corrected
IF (MOREOBS .or. ZERO .ne. 0.) THEN
J = J + 1
LN = -1
CALL FESD((t1-ZERO)/100., -0.001, buffer(j),ln)
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
END IF
END IF
C----------------------------------------------------------------------
C intensity info:
IF (NEEDESD) THEN
C _pd_proc_intensity_total or _pd_meas_intensity_total
ESD = 0
IF (YW .gt. 0) ESD = 1./SQRT(YW)
J = J + 1
LN = -1
CALL FESD(YO, ESD, buffer(j),ln)
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
ELSE
C _pd_meas_intensity_counts
J = J + 1
LN = -1
CALL FESD(YO, -1., buffer(j),ln)
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
END IF
C----------------------------------------------------------------------
C _pd_proc_ls_weight
J = J + 1
IF ( BTEST(ICODE,1) ) THEN
LN = -1
CALL FESD(0.0, 0.0, buffer(j),ln)
ELSE
LN = -1
CALL FESD(YW, -yw/100., buffer(j),ln)
END IF
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
C----------------------------------------------------------------------
C _pd_proc_intensity_calc_bkg
J = J + 1
LN = -1
CALL FESD(YB1, -ESD/10., buffer(j),ln)
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
! if (j .ge. 5) then
! write (IUCIF,'(5(2x,a))')
! 1 (buffer(jj)(:LENCH(buffer(jj))),jj=1,j)
! j = 0
! END IF
C for now ignore fixed background points
! IF (FIXBACK) WRITE(IUCIF,'(6x,A)') '_pd_proc_intensity_fix_bkg'
C----------------------------------------------------------------------
C _pd_calc_intensity_total
J = J + 1
IF(BTEST(ICODE,1)) THEN
buffer(j) = ' .' ! undefined Y(calc)
ELSE
LN = -1
CALL FESD(YC, -esd/10., buffer(j),ln)
END IF
BUFLEN(J) = MAX(BUFLEN(J),LENCH(BUFFER(J)))
C----------------------------------------------------------------------
C write the line to the scratch file
write (IUSCRT,'(9A)') (buffer(jj),jj=1,j)
JMAX = J
END IF
END DO
REWIND(IUSCRT)
C copy from the scratch file to the table
DO I=1,NPOINT
READ(IUSCRT,'(9A)') (BUFFER(JJ),JJ=1,JMAX)
TEXT = ' '
LN = 1
DO JJ=1,JMAX
IF (LN + BUFLEN(JJ) .GT. 80) THEN
WRITE(IUCIF,'(A)') TEXT(:LENCH(TEXT))
TEXT = ' '
LN = 5
END IF
TEXT(LN+1:) = BUFFER(JJ)(1:BUFLEN(JJ))
LN = LN + BUFLEN(JJ) + 2
END DO
WRITE(IUCIF,'(A)') TEXT(:LENCH(TEXT))
END DO
CLOSE(IUSCRT,STATUS='DELETE')
IF (.not. MOREOBS) THEN
write (text,'(I9)') npoint
CALL WRVAL(IUCIF, '_pd_meas_number_of_points', text)
END IF
write (text,'(I9)') npoint
CALL WRVAL(IUCIF, '_pd_proc_number_of_points', text)
CLOSE(IUPRF)
RETURN
END