PROGRAM SMOOTH ! gfortran -O3 -o smoothsod smoothsod.f90 ! last revisions 08jan89 JSK, 04oct93, 24feb94 UFK ! ************************************************************ ! SMOOTH POLE FIGURES or OD sections BY GAUSSIAN DIGITAL FILTER ! ADJUSTED FOR COMPRESSION IN THE CENTER. ! Written by JOHN KALLEND, OCT 1987, Glenwood IL ! (C) Copyright 1987 John S. Kallend ! Permission is hereby granted to copy, use and modify ! this program. Please feel free to give it to your friends! ! ************************************************************ ! IMPLICIT INTEGER*2 (I-N) LOGICAL CPF , not_end 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./ ! DATA not_end/.true./ logical l_com_line character arg_line*20 bs=char(8) CPF = .FALSE. not_end = .true. l_com_line = .false. NREF=1 CALL GetArg(1,fname) ! print*,'debug: input filename = ',fname if(fname.ne.'') l_com_line = .true. ! detect an argument for the input file name; look for user input if not present CALL GetArg(2,arg_line) ! print*,'debug: arg_line = ',arg_line read(arg_line,"(g10.3)") wdth ! 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 end if 1011 FORMAT(A) 1012 FORMAT(A) j=40 do 106 i=40,1,-1 ! write(*,*) fname(i:i) if(fname(i:i).eq.'.') then j=i ! goto 107 exit end if 106 end do 107 continue ! write(*,*) 'j= ',j if(j.eq.1) then ! write(*,*) 'Check, j= ',j spec=fname(1:8) fname=spec//'.cod' else spec=fname(1:j-1) endif ! 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') ! PRINT 1010,' Smoothing range in degs.(w/dec.pt.)--> 2.5',bs,bs,bs 1010 format(/4a) ! 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 ! wdth=7.5 CALL SMARRAY(SMF,WDTH) ! pause ' continue?' print * print 1013,' Making ',fname 1013 FORMAT(2A/) ! **************************** begin loop do while ( not_end ) 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 ! 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 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' ! 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 ! ! 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 do 100 i=1,nrows do 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) ! write(*,*) 'Input: I= ',i ! print 2000,(ICOUNT(i,j),j=ist,ist+ncolr) 2000 format(1x,19i4) end do 100 end do 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 j=1,72 icount(i,j)=icount(i,j)*scale end do 150 end do endif endif WRITE(9,1001) hkl,PINC,PMAX,AINC,AZMAX,IW,JW,IPER,iavg,cstuff do 200 i=1,nrows do 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 ! write(*,*) 'Smoothed: I= ',i ! print 2000,(ICOUNT(i,j),j=ist,ist+ncolr) write(9,2000)(ICOUNT(i,j),j=ist,ist+ncolr) end do 200 end do WRITE(9,*) READ(3,1011,END=35) NREF=NREF+1 GOTO 233 ! ******** end loop END do 35 CLOSE(UNIT=3) print *,' ' ! RETURN call exit ! stand alone program now 55 PRINT*,' CANNOT OPEN FILE ',FNAME stop end program SMOOTH 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 K=NQUADS+1,4 IST=(K-1)*IPQUAD+2 LST=IST+IPQUAD-1 IF(LST.GT.72)LST=72 DO 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) END DO END DO 100 END DO END IF RETURN END SUBROUTINE AZFILL ! ! SUBROUTINE SMARRAY(SMF,WDTH) ! WDTH is the width of the smoothing filter ! SMF is the array that contains the filter values ! IMPLICIT INTEGER*2 (I-N) REAL dif,SMF(19,-36:36) ! DO 100 I=2,19 do 100 i=1,19 ! DEG=(I-1)*0.087266 ! SDEG=SIN(DEG) ! 5*PI/180 ! RWDTH=WDTH/SDEG RWDTH=WDTH ! special version: adr: no adjustment of width for declination angle ! dec 01 ! 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 ! print *,k,smf(i,k) 200 end do TOT=0. ! DO 210 K=-35,36 DO 210 K=-36,36 ! was -35 to +36, WHY???? TOT=TOT+SMF(I,K) 210 end do DO 220 K=-36,36 SMF(I,K)=SMF(I,K)/TOT 220 end do 100 end do ! ! now for special treatment of the first row ! DO 230 K=-36,36 !230 SMF(1,K)=1./72 ! ! write(*,*) 'SMF, smooth array:' DO 300 I=1,19 SMF(I,-36)=SMF(I,-36)/2. SMF(I,36)=SMF(I,-36) ! write(*,"('i= ',i3,9(1x,9f8.3))") i,(SMF(I,K),K=-36,36) ! write(*,"('i= ',i3,1x,9f8.3)") i,(SMF(I,K),K=-35,36) 1000 FORMAT(1X, 9F8.3) 300 end do RETURN END SUBROUTINE SMARRAY ! ! SUBROUTINE SMAZ(SMF,ICOUNT,imax) ! IMPLICIT INTEGER*2 (I-N) INTEGER*2 ICOUNT(19,72) REAL SMF(19,-36:36),SMCT(72) imax=0 ! write(*,*) 'ICOUNT: ' DO 100 I=1,19 ! write(*,*) 'unsmoothed result:' ! 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 end do SMCT(K)=CCNT 150 end do DO 110 K=1,72 ICOUNT(I,K)=SMCT(K)+.5 if(icount(i,k).gt.imax) imax=icount(i,k) 110 end do ! write(*,*) 'Smoothed result:' ! PRINT 1000,(ICOUNT(I,K),K=1,72) 1000 FORMAT(9I8) ! write(*,*) 'imax= ',imax 100 end do RETURN END SUBROUTINE SMAZ ! 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) ! CPF is the flag for complete pole figures 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 END IF IF(INDX.GT.19)THEN INDX=38-INDX KK=K-36 IF(KK.LT.1)KK=KK+72 END IF IF(.NOT.CPF .AND. INDX.GT.17)INDX=17 CCNT=ICOUNT(INDX,KK)*SMF(19,IX)+CCNT 180 end do ncount(i,k)=CCNT 150 end do 100 end do ! do 200 i=1,imax do k=1,72 icount(i,k)=ncount(i,k) ! if(icount(i,k).lt.0)icount(i,k)=0 ! remove line above to allow for negative counts in difference files! end do 200 end do ! ! write(*,*) 'Current frame: ' ! do 950, i=1,19 ! write(*,900) (icount(i,k),k=1,19) !900 format(1x,19i4) !950 continue ! RETURN END SUBROUTINE SMPOL