program scorrel implicit none real*8 x,y,a,c,x2,y2,xy,si,sj,al,fik,wmax,wmin,dw, 1dx,dy integer*4 ns,np,il,is,i,j,k,l,ic,n,kmax,i1,i2,j1,j2 character*80 fn,s80 real*8,allocatable::s(:,:),xs(:),sf(:,:),cp(:,:),av(:), 1craw(:,:) logical lnorm,sav write(6,600) 600 format( 1' Correlation of spectral intensities in a series of spectra',/, 1/, 1 ' Input: SP.LST list of spectra',/ 1 ' SCOR.OPT options',/ 1 ' Output: COR.TXT 2D map') call readopt(wmin,wmax,kmax,dw,lnorm,sav) open(9,file='SP.LST',status='old') ns=0 1 read(9,90,end=99,err=99)fn 90 format(a80) ns=ns+1 if(ns.eq.1)then np=0 is=0 il=0 open(10,file=fn,status='old') write(6,*)fn 2 read(10,90,end=88,err=88)s80 il=il+1 ic=0 read(s80,*,err=77,end=77)x,y ic=ic+1 77 if(ic.ne.0)np=np+1 if(ic.ne.0.and.np.eq.1)is=il goto 2 88 close(10) write(6,601)il,np,is 601 format(' First spectrum:',i6,' lines,',i6, 1 ' points, start on line',i6,'.') endif goto 1 99 close(9) write(6,*)ns,' spectra in SP.LST' c skip some points for safe reading: il=il+1 np=np-2 allocate(s(ns,np),xs(np),av(np)) open(9,file='SP.LST') do 3 i=1,ns read(9,90)fn open(10,file=fn) do 31 j=1,is-1 31 read(10,*) do 32 j=is,is+np-1 32 read(10,*)xs(j-is+1),s(i,j-is+1) 3 close(10) close(9) write(6,*)' Spectra loaded' call dav(av,np,ns,wmin,wmax,xs,s) if(lnorm)call norm(xs,s,ns,np,wmin,wmax) c pre-calculate basis al=wmax-wmin allocate(sf(kmax,np)) do 5 k=1,kmax do 5 i=1,np 5 sf(k,i)=fik(k,xs(i),wmin,al) write(6,*)' basis precalculated' allocate(cp(kmax,kmax),craw(np,np)) cp=0.0d0 craw=0.0d0 a=dble(ns) c calculate Fourier coefficients do 4 i=1,np-1 dx=xs(i+1)-xs(i) do 4 j=1,np-1 dy=xs(j+1)-xs(j) if(xs(i).ge.wmin.and.xs(i).le.wmax)then if(xs(j).ge.wmin.and.xs(j).le.wmax)then x=0.0d0 y=0.0d0 x2=0.0d0 y2=0.0d0 xy=0.0d0 do 41 is=1,ns si=s(is,i) sj=s(is,j) x =x +si y =y +sj x2=x2+si*si y2=y2+sj*sj 41 xy=xy+si*sj c=(a*xy-x*y)/dsqrt((a*x2-x**2)*(a*y2-y**2)) if(sav)c=av(i)*av(j)*c craw(i,j)=c do 42 k=1,kmax do 42 l=1,kmax 42 cp(k,l)=cp(k,l)+c*sf(k,i)*sf(l,j)*dx*dy endif endif 4 continue write(6,*)' coefficients done' c write raw map n=int(al/dw) open(9,file='RAW.TXT') x=wmin-dw do 8 i=1,n x=x+dw call find(i1,i2,x,np,xs) y=wmin-dw do 8 j=1,n y=y+dw call find(j1,j2,y,np,xs) c=0.0d0 if(i1.gt.0.and.i2.gt.0.and.j1.gt.0.and.j2.gt.0)then c=(craw(i1,j1)+craw(i1,j2)+craw(i2,j1)+craw(i2,j2))/4.0d0 endif 8 write(9,98)x,y,c close(9) write(6,*)' RAW.TXT written' c write smooth result open(9,file='COR.TXT') c pre-calculate basis deallocate(sf) allocate(sf(kmax,n)) do 6 k=1,kmax x=wmin-dw do 6 i=1,n x=x+dw 6 sf(k,i)=fik(k,x,wmin,al) write(6,*)' basis II precalculated' x=wmin-dw do 7 i=1,n x=x+dw y=wmin-dw do 7 j=1,n y=y+dw c=0.0d0 do 51 k=1,kmax do 51 l=1,kmax 51 c=c+cp(k,l)*sf(k,i)*sf(l,j) 7 write(9,98)x,y,c 98 format(2F9.2,f10.4) close(9) write(6,*)' COR.TXT written' end subroutine readopt(wmin,wmax,kmax,dw,lnorm,sav) implicit none real*8 dw,wmin,wmax integer*4 kmax character*3 key logical lnorm,sav kmax=100 wmin=200 wmax=3000 dw=30 lnorm=.true. sav=.true. open(9,file='SCOR.OPT') 1 read(9,90,end=99,err=99)key 90 format(a3) if(key.eq.'KMA')read(9,*)kmax if(key.eq.'WMI')read(9,*)wmin if(key.eq.'WMA')read(9,*)wmax if(key.eq.'DEL')read(9,*)dw if(key.eq.'NOR')read(9,*)lnorm if(key.eq.'SAV')read(9,*)sav goto 1 99 close(9) return end c ------------------------------------------------------------ function fik(k,x,xs,al) c k: 1 2 3 4 c cos and sin const sin(2pix/L) cos(2pix/L) sin(3pix/L) c k ...frequency (k=1 ... zero) c x ... variable c al ... interval or HO frequency c xs ... starting x for plane waves (origin) implicit none real*8 fik,pi,x,al,sf,y,xs integer*4 k pi=4.0d0*datan(1.0d0) y=x-xs if(k.eq.1)then sf=1.0d0/dsqrt(al) else sf=dsqrt(2.0d0/al) endif if(mod(k,2).ne.0)then fik=sf*cos(dble(2*((k-1)/2))*pi*y/al) else fik=sf*sin(dble(2*(k/2))*pi*y/al) endif return end c ------------------------------------------------------------ subroutine norm(xs,s,ns,np,wmin,wmax) implicit none integer*4 ns,np,i,j real*8 s(ns,np),wmin,wmax,xs(*),an,x do 1 i=1,ns an=0.0d0 do 11 j=1,np x=xs(j) 11 if(x.ge.wmin.and.x.le.wmax)an=an+dabs(s(i,j)) do 1 j=1,np 1 s(i,j)=s(i,j)/an write(6,*)'Spectra were normalized' return end subroutine find(i1,i2,x,np,xs) implicit none integer*4 i1,i2,np,i real*8 x,xs(*) i1=0 i2=0 do 1 i=1,np-1 if(xs(i).le.x.and.xs(i+1).ge.x)then i1=i i2=i+1 return endif 1 continue return end subroutine dav(av,np,ns,wmin,wmax,xs,s) implicit none integer*4 np,ns,i,j real*8 av(np),s(ns,np),sm,wmin,wmax,xs(*) c average normalized spectrum av=0.0d0 do 1 i=1,ns do 1 j=1,np 1 av(j)=av(j)+s(i,j) sm=0.0d0 do 33 j=1,np if(xs(j).ge.wmin.and.xs(j).le.wmax)then if(dabs(av(j)).gt.sm)sm=dabs(av(j)) endif 33 continue do 34 j=1,np 34 av(j)=av(j)/sm open(9,file='NORM.PRN') do 2 i=1,np 2 write(9,90)xs(i),av(i) 90 format(f12.2,f12.6) close(9) write(6,*)'NORM.PRN average signal' return end