PROGRAM SMOOTH C last revisions 08jan89 JSK, 04oct93, 24feb94 UFK 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 ************************************************************ C IMPLICIT INTEGER*2 (I-N) LOGICAL*1 CPF CHARACTER FNAME*40,DAT*9,title*43,IPFC*1,cstuff*15,hkl*5,spec*20 character line*80,bs INTEGER*2 ICOUNT(19,72),H,K,L,IPER(3),j REAL SMF(19,-36:36) DATA CPF/.FALSE./ logical l_com_line character arg_line*20 bs=char(8) C NREF=1 l_com_line = .false. CALL GetArg(1,fname) c print*,'debug: input filename = ',fname if(fname.ne.'') l_com_line = .true. c detect an argument for the input file name; look for user input if not present CALL GetArg(2,arg_line) c print*,'debug: arg_line = ',arg_line read(arg_line,"(g10.3)") wdth c print*,'debug: wdth = ',wdth PRINT *,' Gaussian smoothing of SOD/COD data' PRINT *,' Program by John Kallend (C) 1987' PRINT *,' Edited 2001 by ADRollett' print *,' to remove dependence of width on decl. angle' print *,' and to add command line arguments, 2004' print * if(.not.l_com_line) then 23 PRINT *,'Input file (with .ext, default .cod): ' READ 1011,FNAME endif 1011 FORMAT(A) 1012 FORMAT(A) j=40 do 106 i=40,1,-1 c write(*,*) fname(i:i) if(fname(i:i).eq.'.') then j=i goto 107 endif 106 continue 107 continue c write(*,*) 'j= ',j if(j.eq.1) then c write(*,*) 'Check, j= ',j spec=fname(1:8) fname=spec//'.cod' else spec=fname(1:j-1) endif c write(*,*) 'spec= ',spec OPEN(FILE=FNAME,UNIT=3,STATUS='OLD',ERR=55) if(fname(j+1:j+1).eq.'S'.or.fname(j+1:j+1).eq.'s') then fname=spec(1:j-1)//'.smh' elseif(fname(j+1:j+1).eq.'C'.or.fname(j+1:j+1).eq.'c') then fname=spec(1:j-1)//'.cmh' else fname=spec(1:j-1)//'.gmh' endif OPEN(UNIT=9,FILE=fname,status='unknown') c PRINT 1010,' Smoothing range in degs.(w/dec.pt.)--> 2.5',bs,bs,bs 1010 format(/4a) c READ '(f4.1)',WDTH if(.not.l_com_line) then write(*,*) 'Enter a smoothing range in degrees [7.5]' READ(*,*) WDTH endif if(wdth.eq.0.) wdth=2.5 c wdth=7.5 CALL SMARRAY(SMF,WDTH) c pause ' continue?' 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) c if(nref.eq.1) THEN c print * c print *,line c IF(PMAX.EQ.90.) THEN c PRINT '(/2a)', c #' Will treat as INCOMPLETE pole figures to 80 deg., OK? Y',bs c READ 1011,IPFC c IF(IPFC.EQ.'N' .OR. IPFC.EQ.'n')CPF=.TRUE. c ENDIF c ENDIF CPF=.TRUE. 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 C ALL DATA TO 90 DEGREES C ncolmn=azmax/AINC/nquads+.01 maxaz=nquads*ncolmn if(nquads.ne.4)maxaz=maxaz+1 maxpo=PMAX/PINC+.1+iw 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 write(*,*) 'Input: I= ',i c print 2000,(ICOUNT(i,j),j=ist,ist+ncolr) 2000 format(1x,19i4) 100 continue if(azmax.lt.360.)CALL AZFILL(IW,NQUADS,MAXPO,MAXAZ,ICOUNT) CALL SMPOL(CPF,SMF,ICOUNT,MAXPO) CALL SMAZ(SMF,ICOUNT,imax) if(iavg.lt.100) then scale=100./float(iavg) amax=imax*scale if(amax.le.9999.) then iavg=100 do 150 i=1,19 do 150 j=1,72 icount(i,j)=icount(i,j)*scale 150 continue endif endif 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 write(*,*) 'Smoothed: I= ',i 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) print *,' ' c RETURN call exit ! stand alone program now 55 PRINT*,' CANNOT OPEN FILE ',FNAME stop END SUBROUTINE AZFILL(IW,NQUADS,MAXPO,MAXAZ,ICOUNT) IMPLICIT INTEGER*2 (I-N) INTEGER*2 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 C SUBROUTINE SMARRAY(SMF,WDTH) c WDTH is the width of the smoothing filter c SMF is the array that contains the filter values IMPLICIT INTEGER*2 (I-N) REAL dif,SMF(19,-36:36) c DO 100 I=2,19 do 100 i=1,19 c DEG=(I-1)*0.087266 c SDEG=SIN(DEG) C 5*PI/180 c RWDTH=WDTH/SDEG RWDTH=WDTH c special version: adr: no adjustment of width for declination angle c dec 01 c PRINT *,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) ELSE SMF(I,K)=EXP(-DIF*DIF) SMF(I,-K)=SMF(I,K) ENDIF c print *,k,smf(i,k) 200 continue TOT=0. c DO 210 K=-35,36 DO 210 K=-36,36 c was -35 to +36, WHY???? TOT=TOT+SMF(I,K) 210 CONTINUE DO 220 K=-36,36 SMF(I,K)=SMF(I,K)/TOT 220 CONTINUE 100 CONTINUE c c now for special treatment of the first row c DO 230 K=-36,36 c230 SMF(1,K)=1./72 c c write(*,*) 'SMF, smooth array:' DO 300 I=1,19 SMF(I,-36)=SMF(I,-36)/2. SMF(I,36)=SMF(I,-36) c write(*,"('i= ',i3,9(1x,9f8.3))") i,(SMF(I,K),K=-36,36) c write(*,"('i= ',i3,1x,9f8.3)") i,(SMF(I,K),K=-35,36) 1000 FORMAT(1X, 9F8.3) 300 CONTINUE RETURN END C C SUBROUTINE SMAZ(SMF,ICOUNT,imax) IMPLICIT INTEGER*2 (I-N) INTEGER*2 ICOUNT(19,72) REAL SMF(19,-36:36),SMCT(72) imax=0 c write(*,*) 'ICOUNT: ' DO 100 I=1,19 c write(*,*) 'unsmoothed result:' c write(*,"(18i4)") (ICOUNT(I,K),K=1,72) 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+ICOUNT(I,INDX)*SMF(I,IX) endif 160 CONTINUE SMCT(K)=CCNT 150 CONTINUE DO 110 K=1,72 ICOUNT(I,K)=SMCT(K)+.5 if(icount(i,k).gt.imax) imax=icount(i,k) 110 CONTINUE c write(*,*) 'Smoothed result:' c PRINT 1000,(ICOUNT(I,K),K=1,72) 1000 FORMAT(9I8) c write(*,*) 'imax= ',imax 100 CONTINUE RETURN END C SUBROUTINE SMPOL(CPF,SMF,ICOUNT,MAXPO) IMPLICIT INTEGER*2 (I-N) LOGICAL*1 CPF INTEGER*2 ICOUNT(19,72),H,K,L,IPER(3),ncount(19,72) REAL SMF(19,-36:36) c CPF is the flag for complete pole figures IF(CPF) THEN IMAX=19 ELSE IMAX=MAXPO IF(MAXPO.EQ.19)MAXPO=17 ENDIF c 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 c do 200 i=1,imax do 200 k=1,72 icount(i,k)=ncount(i,k) c if(icount(i,k).lt.0)icount(i,k)=0 c remove line above to allow for negative counts in difference files! 200 continue c c write(*,*) 'Current frame: ' c do 950, i=1,19 c write(*,900) (icount(i,k),k=1,19) c900 format(1x,19i4) c950 continue c RETURN END