PROGRAM SMOOTH C last revisions 08jan89 JSK, 04oct93, 24feb94 UFK c Sept 06, ADR, adapt to unix C ************************************************************ C SMOOTH POLE FIGURES or OD sections BY GAUSSIAN DIGITAL FILTER C ADJUSTED FOR COMPRESSION IN THE CENTER. C Written by JOHN KALLEND, OCT 1987, Glenwood IL C (C) Copyright 1987 John S. Kallend C Permission is hereby granted to copy, use and modify C this program. Please feel free to give it to your friends! C ************************************************************ IMPLICIT integer (I-N) LOGICAL*1 CPF c CHARACTER FNAME*12,DAT*9,title*43,IPFC*1,cstuff*15,hkl*5,spec*8 CHARACTER DAT*9,title*43,IPFC*1,cstuff*15,hkl*5 character fname*50,spec*48 character line*80,bs integer ICOUNT(19,72),H,K,L,IPER(3) REAL SMF(19,-36:36) DATA CPF/.FALSE./ c CODE:: bs=char(8) NREF=1 PRINT 1011,'Gaussian smoothing of pole figure data' PRINT 1011,' Program by John Kallend (C) 1987' print*,'modified by ADR for unix' print*,'[/path/]smooth data width[2.5] incomplete[Y/N]' if ( iargc() .gt. 0 ) then call getarg ( 1 , fname ) else 23 PRINT 1012,'Input file (with .ext, default .EPF): ' READ 1011,FNAME 1011 FORMAT(A) c 1012 FORMAT(A\) 1012 FORMAT(A) endif j = 0 c$$$ do 106 i=1,9 c$$$ if(fname(i:i).eq.'.')j=i c$$$ 106 continue j = index(fname,'.') if(j.eq.0) j = index(fname,' ') print*,'dot position at ',j c if(j.eq.10) then if(j.ge.46) then spec=fname(1:46) fname=spec//'.EPF' else spec=fname(1:j-1) endif OPEN(FILE=FNAME,UNIT=3,STATUS='OLD',ERR=55) fname=spec(1:j-1)//'.MPF' OPEN(UNIT=9,FILE=fname) if ( iargc() .gt. 1 ) then call getarg ( 2 , line ) read(line,*) wdth else PRINT"(' Smoothing range in degs.(w/dec.pt.)--> 2.5',a,a,a,$)" 1 ,bs,bs,bs c 1010 format(/4a\) 1010 format(/4a) READ '(f4.1)',WDTH endif if(wdth.eq.0) wdth=2.5 CALL SMARRAY(SMF,WDTH) print * print 1013,' Making ',fname 1013 FORMAT(2A/) c **************************** begin loop 233 READ (3,1011,end=35)line title=line WRITE(9,1021)title//', Smoothed',WDTH,' deg.'//line(64:73) 1021 format(a,f5.1,a) READ (3,1001) hkl,PINC,PMAX,AINC,AZMAX,IW,JW,IPER,iavg,cstuff 1001 FORMAT(a5,4F5.1,5I2,i5,a15) if(nref.eq.1) THEN print * print *,line IF(PMAX.EQ.90.) THEN c PRINT '(/2a\)', if ( iargc() .gt. 2 ) then call getarg ( 3 , ipfc ) else PRINT '(/2a,$)', &' Will treat as INCOMPLETE pole figures to 80 deg., OK? Y',bs READ 1011,IPFC IF(IPFC.EQ.'N' .OR. IPFC.EQ.'n') CPF = .TRUE. endif ENDIF ENDIF print 1000, hkl,PINC,PMAX,AINC,AZMAX,IW,JW,IPER,iavg,cstuff 1000 FORMAT(1x,a5,4F5.1,5I2,i5,a15) IF(JW*iw.eq.0) stop'cannot handle iw=0 or jw=0 data format' C IF(PMAX.LT.80.)STOP 'INSUFFICIENT DATA - need 80 degrees' IF(pinc.ne.5.0) stop'can only handle increments of 5 deg' nquads=azmax/90.+.01 nrows=90./PINC+.01+iw C ALL DATA TO 90 DEGREES ncolmn=azmax/AINC/nquads+.01 maxaz=nquads*ncolmn if(nquads.ne.4) maxaz = maxaz+1 maxpo=PMAX/PINC+.1+iw c print*,'reading in data ' do 100 i=1,nrows do 100 k=1,nquads NCOLR=ncolmn-1 if(k.ne.4 .and. iw.eq.1 .and. k.eq. nquads)ncolr=ncolmn ist=(k-1)*ncolmn+1 read (3,2000)(ICOUNT(i,j),j=ist,ist+ncolr) c print 2000,(ICOUNT(i,j),j=ist,ist+ncolr) 2000 format(1x,19i4) 100 continue if(azmax.lt.360.) then CALL AZFILL(IW,NQUADS,MAXPO,MAXAZ,ICOUNT) c print*,'ran AZFILL' endif CALL SMPOL(CPF,SMF,ICOUNT,MAXPO) CALL SMAZ(SMF,ICOUNT,imax) if(iavg.lt.100) then scale=100./float(iavg) amax=imax*scale print*,' scale, amax = ',scale,amax if(amax.le.9999.) then iavg=100 do 150 i=1,19 do 150 j=1,72 icount(i,j) = int(float(icount(i,j))*scale) 150 continue else stop 'AMAX too large?!' endif else print*,'IAVG >= 100' endif c print*,'scaled data - writing out to MPF' WRITE(9,1001) hkl,PINC,PMAX,AINC,AZMAX,IW,JW,IPER,iavg,cstuff do 200 i=1,nrows do 200 k=1,nquads NCOLR=ncolmn-1 if(k.ne.4 .and. iw.eq.1 .and. k.eq. nquads)ncolr=ncolmn ist=(k-1)*ncolmn+1 c print 2000,(ICOUNT(i,j),j=ist,ist+ncolr) write(9,2000)(ICOUNT(i,j),j=ist,ist+ncolr) 200 continue WRITE(9,*) READ(3,1011,END=35) NREF=NREF+1 GOTO 233 c ******** end loop 35 CLOSE(UNIT=3) c print *,' ' print *,' done!' RETURN 55 PRINT 1013,' CANNOT OPEN FILE ',FNAME call exit(1) ! PRINT 1011 ! GOTO 23 END SUBROUTINE AZFILL(IW,NQUADS,MAXPO,MAXAZ,ICOUNT) IMPLICIT integer (I-N) integer ICOUNT(19,72) IPQUAD=(MAXAZ-1)/NQUADS+.01 IF(IW.EQ.1)THEN DO 100 I=1,MAXPO DO 100 K=NQUADS+1,4 IST=(K-1)*IPQUAD+2 LST=IST+IPQUAD-1 IF(LST.GT.72)LST=72 DO 100 J=IST,LST IF(K.EQ.2)JAC=IPQUAD*2+2-J IF(K.EQ.3 .OR. K.EQ.4)JAC=J-2*IPQUAD ICOUNT(I,J)=ICOUNT(I,JAC) 100 CONTINUE ENDIF RETURN END C SUBROUTINE SMARRAY(SMF,WDTH) IMPLICIT integer (I-N) REAL SMF(19,-36:36) c print*,'starting SMARRAY, WDTH= ',WDTH DO 100 I = 2,19 DEG=(I-1)*0.087266 SDEG=SIN(DEG) C 5*PI/180 RWDTH = WDTH/SDEG c PRINT *,'i,rwdth,sdeg= ',I,RWDTH,SDEG DO 200 K=0,36 DIF = float(K)*5./RWDTH IF(DIF.GT.8) THEN SMF(I,K) = 0. SMF(I,-K)=SMF(I,K) c JSK did not include this line but we cannot assume zero initialization! ELSE SMF(I,K)=EXP(-DIF*DIF) SMF(I,-K)=SMF(I,K) ENDIF c print *,k,smf(i,k) 200 continue c print *,'SMF: ',(smf(i,k),k=0,36) TOT = 0. DO 210 K=-35,36 TOT=TOT+SMF(I,K) 210 CONTINUE c print*,'I, TOT = ',I,TOT c if(tot.eq.0.) stop 'TOT = 0!' DO 220 K=-36,36 SMF(I,K) = SMF(I,K)/TOT 220 CONTINUE 100 CONTINUE DO 230 K = -36,36 230 SMF(1,K) = 1./72 DO 300 I=1,19 c print*,'i= ',i SMF(I,-36)=SMF(I,-36)/2. SMF(I,36)=SMF(I,-36) c PRINT 1000,(SMF(I,K),K=-35,36) 1000 FORMAT(1X, 9F8.3) 300 CONTINUE RETURN END C SUBROUTINE SMAZ(SMF,ICOUNT,imax) IMPLICIT integer (I-N) integer ICOUNT(19,72) REAL SMF(19,-36:36),SMCT(72) c CODE:: c print*,'starting SMAZ' c PRINT 1000,(ICOUNT(I,K),K=1,72) imax=0 DO 100 I=1,19 DO 150 K=1,72 CCNT=0. DO 160 IX=-36,36 if(smf(i,ix).gt.0.005)then INDX=K+IX IF(INDX.LT.1)INDX=INDX+72 IF(INDX.GT.72)INDX=INDX-72 CCNT=CCNT+float(ICOUNT(I,INDX))*SMF(I,IX) endif 160 CONTINUE SMCT(K)=CCNT 150 CONTINUE DO 110 K=1,72 ICOUNT(I,K) = int(SMCT(K)+.5) if(icount(i,k).gt.imax) imax = icount(i,k) 110 CONTINUE c PRINT 1000,(ICOUNT(I,K),K=1,72) 1000 FORMAT(18I4) 100 CONTINUE RETURN END C SUBROUTINE SMPOL(CPF,SMF,ICOUNT,MAXPO) IMPLICIT integer (I-N) LOGICAL*1 CPF integer ICOUNT(19,72),H,K,L,IPER(3),ncount(19,72) REAL SMF(19,-36:36) ! print*,'starting SMPOL' ! PRINT 1000,(ICOUNT(I,K),K=1,72) 1000 FORMAT(18I4) IF(CPF) THEN IMAX=19 ELSE IMAX=MAXPO IF(MAXPO.EQ.19)MAXPO=17 ENDIF DO 100 K=1,72 DO 150 I=1,IMAX CCNT=0. DO 180 IX=-6,6 KK=K INDX=I+IX IF(INDX.LT.1)THEN INDX=2-INDX KK=K-36 IF(KK.LT.1)KK=KK+72 ENDIF IF(INDX.GT.19)THEN INDX=38-INDX KK=K-36 IF(KK.LT.1)KK=KK+72 ENDIF IF(.NOT.CPF .AND. INDX.GT.17)INDX=17 180 CCNT=ICOUNT(INDX,KK)*SMF(19,IX)+CCNT ncount(i,k)=CCNT 150 CONTINUE 100 CONTINUE do 200 i=1,imax do 200 k=1,72 icount(i,k)=ncount(i,k) if(icount(i,k).lt.0)icount(i,k)=0 200 continue RETURN END