program specomp implicit none integer*4 ns0,np0 parameter (ns0=2000,np0=9000) character*80 sn(ns0), s80 character*10 s10 integer*4 is,i,j,js,IERR,ii,ndim real*8 as(ns0),ks(ns0),tem,bs(ns0),norm1,normj,so(ns0,ns0), 1vo(ns0),e(ns0,2*ns0),ai(ns0,ns0),TOL,coe(ns0),sum,ANG(3),FU, 2cavused,st(np0),ps,pss,coeold(ns0),err,serr integer n1,n2,ITER,nall(ns0),ns,jsc real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0),RTOL,RMS,a,k,b, 1sp(ns0),sallx(ns0,np0),sally(ns0,np0),sfit(ns0,np0),rs COMMON/MATVAR/s1x,s1y,s2x,s2y,n1,n2,RTOL,ITER,RMS,rs,a,k,b, 1sallx,sally,nall,ns,sp,sfit,jsc real*8 wmin,wmax,aa,kk,bb,al,kin,bin,cav,xs0,xs1,xs2 integer*4 npar,iwr,icol,iinv,ino,normc,ertype,iraw,adapt,xscale, 1ienergy common/opt/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino,normc,ertype,iraw,adapt,cav,xscale,xs0,xs1,xs2,ienergy c write(6,6000) 6000 format(/,' Specomp compares spectra listed in SP.LST',/, 1 ' to the first spectrum in the list.',/) call iopar c open(9,file='SP.LST') ns=0 1 read(9,9000,end=99,err=99)s80 9000 format(a80) if(s80(1:1).eq.'#')goto 1 ns=ns+1 if(ns.gt.ns0)call report('too many spectra') sn(ns)=s80 goto 1 99 close(9) if(ertype.eq.0)then write(6,*)' abs(A-B) used as error' else if(ertype.eq.1)then write(6,*)' sqrt((A-B)**2) used as error' else if(ertype.eq.2)then write(6,*)' |A-B|/(|B|) used as error' else if(ertype.eq.3)then write(6,*)' (A-B)**2/(|A||B|) used as error' else if(ertype.eq.4)then write(6,*)' A.B/(|A||B|) used as error' else write(6,*)ertype call report('Unknown error option') endif endif endif endif endif if(icol.gt.0)write(6,*)'Collective spectra fitting' if(icol.le.0)write(6,*)'Individual spectra fitting' if(normc.eq.0)then write(6,*)'Normalize all spectra' if(ino.eq.0)write(6,*)'Square normalizing.' if(ino.ne.0)write(6,*)'Absolute intensity normalizing.' else write(6,*)'Raw input spectra used' norm1=1.0d0 normj=1.0d0 endif write(6,6003)ns 6003 format(i4,' spectra',/,/, 1'spec a k b RMS rs ITER') c c read decomposed spectrum: is=1 call ls(s1x,s1y,n1,np0,sn(is)) if(normc.eq.0)call norm(n1,s1x,s1y,norm1,wmin,wmax,ino) c if(icol.gt.0)then call sav(1,n1,s1x,s1y,sallx,sally,nall,ns0,np0) do 7 js=2,ns call ls(s2x,s2y,n2,np0,sn(js)) if(xscale.eq.1)call ss(s2x,n2,xs0,xs1,xs2) if(normc.eq.0)call norm(n2,s2x,s2y,normj,wmin,wmax,ino) 7 call sav(js,n2,s2x,s2y,sallx,sally,nall,ns0,np0) if(adapt.ne.0)then call DOMATRIX else a=1.0d0 k=1.0d0 b=0.0d0 endif as(1)=1.0d0 ks(1)=1.0d0 bs(1)=0.0d0 do 9 js=2,ns as(js)=a ks(js)=k bs(js)=b 9 if(npar.eq.2)bs(js)=0.0d0 else do 2 js=1,ns jsc=js call ls(s2x,s2y,n2,np0,sn(js)) if(xscale.eq.1)call ss(s2x,n2,xs0,xs1,xs2) if(normc.eq.0)call norm(n2,s2x,s2y,normj,wmin,wmax,ino) if(adapt.ne.0)then call DOMATRIX else a=1.0d0 k=1.0d0 b=0.0d0 ANG(1)=a ANG(2)=k ANG(3)=b ITER=0 RMS=FU(ANG) rs=RMS endif write(6,6002)js,a,k,b,RMS,rs,ITER 6002 format(i4,3f8.3,2g15.6,i5) as(js)=a ks(js)=k bs(js)=b if(npar.eq.2)bs(js)=0.0d0 c c write intermediate spectra: if(iwr.gt.0)then write(s10,1010)js 1010 format(i10) do 5 is=1,10 5 if(s10(is:is).ne.' ')goto 6 6 call WRFU(a,k,b,s1x,s2x,s2y,n1,n2,'_'//s10(is:10)//'.prn',norm1) if(iwr.gt.1)call WRS(s2x,s2y,n2,'_'//s10(is:10)//'N.prn') endif c c write(6,*)js,RMS 2 sp(js)=RMS endif c open(9,file='SP.OUT') do 44 i=1,ns 44 write(9,9001)i,sp(i),as(i),ks(i),bs(i),sn(i)(1:30) write(9,*) do 3 i=1,ns do 3 j=i+1,ns if(sp(i).gt.sp(j))then tem=as(i) as(i)=as(j) as(j)=tem tem=bs(i) bs(i)=bs(j) bs(j)=tem tem=ks(i) ks(i)=ks(j) ks(j)=tem tem=sp(i) sp(i)=sp(j) sp(j)=tem s80=sn(i) sn(i)=sn(j) sn(j)=s80 endif 3 continue write(6,6001) 6001 format('SP.OUT written, best spectra:',/, 1' # RMS a k b') do 4 i=1,ns if(i.lt.4)write(6,9001)i,sp(i),as(i),ks(i),bs(i),sn(i)(1:30) 4 write(9,9001)i,sp(i),as(i),ks(i),bs(i),sn(i)(1:30) 9001 format(i4,g12.4,3f8.3,' ',a30) close(9) if(iinv.gt.0)then c if(cav.eq.0.0d0)then cavused=1.0d0/dble(ns-1) else cavused=cav endif write(6,609)al,cavused 609 format(/,' Working alpha: ',e10.4,' cave: ',e10.4,/) call vm(so,ns0,ns) do 19 i=2,ns do 19 j=2,ns if(i.eq.j)so(i-1,j-1)=al do 19 ii=1,n1 19 if(s1x(ii).ge.wmin.and.s1x(ii).le.wmax) 1 so(i-1,j-1)=so(i-1,j-1)+sfit(i,ii)*sfit(j,ii) do 10 i=2,ns vo(i-1)=al*cavused do 10 ii=1,n1 10 if(s1x(ii).ge.wmin.and.s1x(ii).le.wmax) 1 vo(i-1)=vo(i-1)+s1y(ii)*sfit(i,ii) write(6,6011)iraw 6011 format(' iraw = ',i1,', S = sum_i c_i S_i') TOL=1.0d-10 if(iraw.gt.0)then write(6,6009) 6009 format(' raw decomposition') ndim=ns-1 else c decompose first spectrum into the other, sum of coefficients c is one and they are slightly force to 1/N: write(6,6010) 6010 format(' normalized decomposition sum_i c_i = 1') do 191 i=2,ns so(ns,i-1)=1.0d0 191 so(i-1,ns)=1.0d0 vo(ns)=1 ndim=ns endif call INV(ns0,so,ai,ndim,TOL,e,IERR) if(IERR.ne.0)call report('unsuccessful inversion') write(6,*)' Spectra decomposition coefficients:' sum=0.0d0 do 11 i=2,ns coe(i)=0.0d0 do 12 ii=2,ndim+1 12 coe(i)=coe(i)+ai(i-1,ii-1)*vo(ii-1) sum=sum+coe(i) 11 write(6,678)i,coe(i) 678 format(i4,f19.4) if(dabs(sum-1.0d0).gt.1.0d-3)then write(6,*)' Coefficients after renormalization:' do 20 i=2,ns 20 write(6,678)i,coe(i)/sum endif if(ienergy.gt.0)then write(6,*)'Energy-based iterative fit' ndim=ns-1 c initial guess: do 21 i=2,ns coe(i)=1.0d0/dble(ns-1) 21 coeold(i)=coe(i) iter=0 999 iter=iter+1 do 22 ii=1,n1 st(ii)=0.0d0 do 22 i=2,ns 22 st(ii)=st(ii)+coe(i)*sfit(i,ii) serr=0.0d0 do 35 ii=1,n1 35 if(s1x(ii).ge.wmin.and.s1x(ii).le.wmax) 1 serr=serr+dabs(st(ii)-s1y(ii)) write(6,*)'initial spectral error',serr ps=0.0d0 pss=0.0d0 do 23 ii=1,n1 if(s1x(ii).ge.wmin.and.s1x(ii).le.wmax)then ps =ps +st(ii)*st(ii) pss=pss+st(ii)*s1y(ii) endif 23 continue c new matrix: call vm(so,ns0,ns) norm1=pss/ps do 24 i=2,ns do 24 j=2,ns do 25 ii=1,n1 25 if(s1x(ii).ge.wmin.and.s1x(ii).le.wmax) 1 so(i-1,j-1)=so(i-1,j-1)+sfit(i,ii)*sfit(j,ii) 24 so(i-1,j-1)=so(i-1,j-1)*norm1 call INV(ns0,so,ai,ndim,TOL,e,IERR) if(IERR.ne.0)call report('unsuccessful inversion') if(iter.gt.100)call report('too many iterations') sum=0 do 27 i=2,ns coe(i)=0.0d0 do 28 ii=2,ndim+1 28 coe(i)=coe(i)+ai(i-1,ii-1)*vo(ii-1) coe(i)=max(0.0d0,coe(i)) 27 sum=sum+coe(i) write(6,*)' Spectra decomposition coefficients:' err=0.0d0 do 33 i=2,ns coe(i)=coe(i)/sum write(6,678)i,coe(i) err=err+dabs(coe(i)-coeold(i)) 33 coeold(i)=coe(i) serr=0.0d0 do 34 ii=1,n1 34 if(s1x(ii).ge.wmin.and.s1x(ii).le.wmax) 1 serr=serr+dabs(st(ii)-s1y(ii)) write(6,800)iter,err,serr 800 format(' Iteration ',i4,' error: ',g12.4,', spectral error: ', 1 g12.4) if(err.lt.1.0d-10)then write(6,*)'Converged after ',iter,' iterations' else goto 999 endif endif c scale the original spectra and reconstruct the sum: do 13 ii=1,n1 13 s2y(ii)=0.0d0 do 61 i=2,ns do 14 ii=1,n1 s1y(ii)=coe(i)*sfit(i,ii)*norm1 14 s2y(ii)=s2y(ii)+coe(i)*sfit(i,ii) c write intermediate spectra: write(s10,1010)i do 51 is=1,10 51 if(s10(is:is).ne.' ')goto 61 61 if(iwr.gt.1)call WRS(s1x,s1y,n1,'_'//s10(is:10)//'P.prn') open(61,file='COEF.TXT') do 32 i=2,ns 32 write(61,678)i,coe(i) close(61) c open(67,file='fit.prn') do 15 i=1,n1 15 if(s1x(i).ge.wmin.and.s1x(i).le.wmax) 1 write(67,900)s1x(i),s2y(i)*norm1 900 format(f9.2,g13.5) close(67) endif stop end subroutine ss(xs,np,xs0,xs1,xs2) c scale x-scale x' = xs0 + xs1 x^1 + xs2 x^2 implicit none integer*4 i,np real*8 xs(*),xs0,xs1,xs2 do 1 i=1,np 1 xs(i)=xs0+xs(i)*(xs1+xs2*xs(i)) return end subroutine ls(xs,ys,np,np0,name) c loads in the spectrum implicit none integer*4 np0,np real*8 xs(*),ys(*),x,y character*(*) name open(9,file=name) np=0 1 read(9,*,end=99,err=99)x,y np=np+1 if(np.gt.np0)call report('too many points') xs(np)=x ys(np)=y goto 1 99 close(9) return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine norm(n,sx,sy,an,wmin,wmax,ino) integer*4 n,i,ino real*8 an,sx(*),sy(*),wmin,wmax an=0.0d0 if(ino.eq.0)then do 1 i=1,n-1 1 if(sx(i).ge.wmin.and.sx(i).le.wmax) 1 an=an+sy(i)**2*abs(sx(i+1)-sx(i)) an=dsqrt(an) else do 11 i=1,n-1 11 if(sx(i).ge.wmin.and.sx(i).le.wmax) 1 an=an+abs(sy(i))*abs(sx(i+1)-sx(i)) endif if(an.gt.1.0d-24)then do 2 i=1,n 2 sy(i)=sy(i)/an endif return end SUBROUTINE DOMATRIX C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT none integer*4 npar0,ITMAX,np0,ns0 parameter (npar0=3,ITMAX=500,np0=9000,ns0=2000) real*8 P(4,npar0),Y(4),PBAR(npar0),PR(npar0),PRR(npar0), 1FU,FTOL,YPR,YPRR,ANG(npar0) parameter (FTOL=1.0d-4) integer*4 ILO,INHI,IHI,I,J integer n1,n2,ITER,nall(ns0),ns,jsc real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0),RTOL,RMS,a,k,b, 1sp(ns0),sallx(ns0,np0),sally(ns0,np0),sfit(ns0,np0),rs COMMON/MATVAR/s1x,s1y,s2x,s2y,n1,n2,RTOL,ITER,RMS,rs,a,k,b, 1sallx,sally,nall,ns,sp,sfit,jsc real*8 wmin,wmax,aa,kk,bb,al,kin,bin,cav,xs0,xs1,xs2 integer*4 npar,iwr,icol,iinv,ino,normc,ertype,iraw,adapt,xscale, 1ienergy common/opt/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino,normc,ertype,iraw,adapt,cav,xscale,xs0,xs1,xs2, 1ienergy c C STARTING VALUES FOR THE ITERATION: c P(,1)..a, P(,2)..k, P(,3)..b if(npar.gt.npar0)call report('too many parameters') P(1,1)=1.0d0 P(1,2)=1.0d0 P(1,3)=0.0d0 P(2,1)=1.0d0 P(2,2)=1.1d0 P(2,3)=1.0d1 P(3,1)=1.1d0 P(3,2)=1.0d0 P(3,3)=5.0d0 P(4,1)=1.1d0 P(4,2)=1.1d0 P(4,3)=2.0d0 do 1 i=1,4 ANG(1)=P(i,1) ANG(2)=P(i,2) ANG(3)=P(i,3) 1 Y(i)=FU(ANG) c ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-10)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(6,*)' Iteration has not converged !' RETURN ENDIF ITER=ITER+1 DO 55 I=1,npar 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,npar 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,npar PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,npar 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,npar 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,npar 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,npar 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,npar 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,npar 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,npar PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,npar 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END subroutine WRFU(a,k,b,s1x,s2x,s2y,n1,n2,fn,norm1) implicit none INTEGER*4 i,n1,jc,n2 REAL*8 a,k,b,s2x(*),s1x(*),w,kw,s2y(*),norm1 character*(*) fn c real*8 wmin,wmax,aa,kk,bb,al,kin,bin,cav,xs0,xs1,xs2 integer*4 npar,iwr,icol,iinv,ino,normc,ertype,iraw,adapt,xscale, 1ienergy common/opt/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino,normc,ertype,iraw,adapt,cav,xscale,xs0,xs1,xs2, 1ienergy c open(67,file=fn) do 1 i=1,n1 w=s1x(i) if(w.ge.wmin.and.w.le.wmax)then kw=w if(npar.gt.1)kw=k*w if(npar.gt.2)kw=kw+b do 2 jc=1,n2-1 2 if((s2x(jc).le.kw.and.s2x(jc+1).ge.kw).or. 1 (s2x(jc).ge.kw.and.s2x(jc+1).le.kw))goto 99 write(67,900)w,0.0d0 goto 1 99 write(67,900)w,a*s2y(jc)*norm1 900 format(f9.2,' ',g13.5) endif 1 continue close(67) write(6,*)fn,' written' RETURN END subroutine WRS(sx,sy,n,fn) implicit none INTEGER*4 i,n REAL*8 sx(*),sy(*) character*(*) fn c open(67,file=fn) do 1 i=1,n 1 write(67,900)sx(i),sy(i) 900 format(f9.2,' ',g12.5) close(67) write(6,*)fn,' written' RETURN END FUNCTION FU(ANG) implicit none INTEGER*4 i,np0,ns0 parameter (np0=9000,ns0=2000) REAL*8 FU,ANG(*),dd,w,kw,anorm1,anorm2,dda integer*4 jc,j integer n1,n2,ITER,nall(ns0),ns,jsc,id real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0),RTOL,RMS,a,k,b, 1sp(ns0),sallx(ns0,np0),sally(ns0,np0),sfit(ns0,np0),rs COMMON/MATVAR/s1x,s1y,s2x,s2y,n1,n2,RTOL,ITER,RMS,rs,a,k,b, 1sallx,sally,nall,ns,sp,sfit,jsc real*8 wmin,wmax,aa,kk,bb,al,kin,bin,cav,xs0,xs1,xs2 integer*4 npar,iwr,icol,iinv,ino,normc,ertype,iraw,adapt,xscale, 1ienergy common/opt/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino,normc,ertype,iraw,adapt,cav,xscale,xs0,xs1,xs2, 1ienergy c anorm1=0.0d0 anorm2=0.0d0 a=ANG(1) k=ANG(2) b=ANG(3) RMS=0.0d0 if(icol.gt.0)then do 8 j=1,ns 8 sp(j)=0.0d0 endif do 1 i=1,n1 w=s1x(i) if(w.ge.wmin.and.w.le.wmax)then kw=w if(npar.gt.1)kw=k*w if(npar.gt.2)kw=kw+b c c difference to all spectra: if(icol.gt.0)then do 2 j=2,ns jc=1 dd=dabs(kw-sallx(j,jc)) 22 jc=jc+1 if(jc.lt.nall(j).and.dabs(kw-sallx(j,jc)).lt.dd)then dd=dabs(kw-sallx(j,jc)) goto 22 endif jc=jc-1 sfit(j,i)=a*sally(j,jc) sp(j)=sp(j)+dabs(sfit(j,i)-s1y(i)) 2 RMS=RMS+dabs(sfit(j,i)-s1y(i)) c c difference to second spectrum: else c c closest point id=1 dda=dabs(s2x(1)-kw) do 24 jc=1,n2 if(dabs(s2x(jc)-kw).lt.dda)then dda=dabs(s2x(jc)-kw) id=jc endif 24 continue jc=id if(dda.gt.max(dabs(s1x(1)-s1x(2)),dabs(s2x(1)-s2x(2))))then sfit(jsc,i)=0.0d0 else sfit(jsc,i)=a*s2y(jc) endif if(ertype.eq.0)then RMS=RMS+dabs(a*s2y(jc)-s1y(i)) else if(ertype.eq.1)then RMS=RMS+(a*s2y(jc)-s1y(i))**2 else if(ertype.eq.2)then RMS=RMS+dabs(a*s2y(jc)-s1y(i)) anorm1=anorm1+dabs(s1y(i)) else if(ertype.eq.3)then RMS=RMS+(a*s2y(jc)-s1y(i))**2 anorm1=anorm1+s1y(i)**2 anorm2=anorm2+(a*s2y(jc))**2 else if(ertype.eq.4)then RMS=RMS +a*s2y(jc)*s1y(i) anorm1=anorm1+ s1y(i)*s1y(i) anorm2=anorm2+a*s2y(jc)*a*s2y(jc) else write(6,*)ertype call report('unknown error type') endif endif endif endif endif endif c endif 1 continue if(ertype.eq.2)RMS=RMS/anorm1 if(ertype.eq.3)RMS=RMS/dsqrt(anorm1*anorm2) if(ertype.eq.4)RMS=RMS/dsqrt(anorm1*anorm2) c c remember only spectral error: rs=RMS c c add also the parametr burden: RMS=RMS+aa*(1.0d0-a)**2 if(npar.gt.1)RMS=RMS+kk*(kin-k)**2 if(npar.gt.2)RMS=RMS+bb*(b-bin)**2 FU=RMS RETURN END subroutine iopar character*10 s10 real*8 wmin,wmax,aa,kk,bb,al,kin,bin,cav,xs0,xs1,xs2 integer*4 npar,iwr,icol,iinv,ino,normc,ertype,iraw,adapt,xscale, 1ienergy common/opt/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin, 1ino,normc,ertype,iraw,adapt,cav,xscale,xs0,xs1,xs2, 1ienergy logical lex ienergy=0 cav=0.0d0 adapt=0 iraw=0 ertype=0 npar=3 wmin=-10000.0d0 wmax=+10000.0d0 kk=1.0d3 aa=1.0d1 bb=4.0d-3 iwr=0 icol=0 al=0.001d0 iinv=0 kin=1.0d0 bin=0.0d0 ino=0 normc=0 xscale=0 xs0=0.0d0 xs1=1.0d0 xs2=0.0d0 inquire(file='SP.OPT',exist=lex) if(lex)then open(9,file='SP.OPT') 1 read(9,90,end=99,err=99)s10 90 format(a10) if(s10(1:4).eq.'adap')read(9,*)adapt if(s10(1:4).eq.'ener')read(9,*)ienergy if(s10(1:5).eq.'normc')read(9,*)normc if(s10(1:5).eq.'ertyp')read(9,*)ertype if(s10(1:4).eq.'alph')read(9,*)al if(s10(1:4).eq.'wmax')read(9,*)wmax if(s10(1:4).eq.'wmin')read(9,*)wmin if(s10(1:2).eq.'kk')read(9,*)kk if(s10(1:2).eq.'aa')read(9,*)aa if(s10(1:2).eq.'bb')read(9,*)bb if(s10(1:3).eq.'bin')read(9,*)bin if(s10(1:3).eq.'kin')read(9,*)kin if(s10(1:4).eq.'iraw')read(9,*)iraw if(s10(1:4).eq.'npar')read(9,*)npar if(s10(1:4).eq.'iinv')read(9,*)iinv if(s10(1:3).eq.'ino')read(9,*)ino if(s10(1:3).eq.'cav')read(9,*)cav if(s10(1:3).eq.'xsc')read(9,*)xscale if(s10(1:3).eq.'xs0')read(9,*)xs0 if(s10(1:3).eq.'xs1')read(9,*)xs1 if(s10(1:3).eq.'xs2')read(9,*)xs2 c c print control: if(s10(1:3).eq.'iwr')read(9,*)iwr c c collective spectra fitting: if(s10(1:4).eq.'icol')read(9,*)icol goto 1 99 close(9) endif write(6,6000)wmin,wmax,aa,kk,bb,al,kin,bin,npar,ienergy, 1ertype,adapt, 1xscale,cav 6000 format('wmin: ',f8.1,' wmax: ',f8.1,/, 1'parameter restrictions:',/, 2' aa (intenzity) : ',f8.1,/, 3' kk (x-scale multiplier): ',f8.1,/, 4' bb (x-scale shift) : ',f8.4,/,/, 6' al (1/N, decomposition): ',f8.4,/,/, 7' ki (syst. x-multiplier): ',f8.1,/, 8' bi (syst. x-shift) : ',f8.1,/, 5' npar : ',i8,/, 5' energy iterative fit : ',i8,/, 5' ertype : ',i8,/, 5' adapt : ',i8,/, 5' xscale : ',i8,/, 5' Cave : ',e8.3,/) if(xscale.eq.1)write(6,6001)xs0,xs1,xs2 6001 format('xs0 : ',f8.3,/, 1 'xs1 : ',f8.3,/, 1 'xs2 : ',f8.2,/) return end subroutine sav(is,ns,sx,sy,sallx,sally,nall,ns0,np0) integer*4 is,i,ns,nall(*),ns0,np0 real*8 sx(*),sy(*),sallx(ns0,np0),sally(ns0,np0) nall(is)=ns do 1 i=1,ns sallx(is,i)=sx(i) 1 sally(is,i)=sy(i) return end SUBROUTINE INV(nr,a,ai,n,TOL,e,IERR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(nr,nr),a(nr,nr),e(nr,2*nr) C 10000 IERR=0 DO 1 ii=1,n DO 1 jj=1,n e(ii,jj)=a(ii,jj) 1 e(ii,jj+n)=0.0D0 do 13 ii=1,n 13 e(ii,ii+n)=1.0D0 c DO 2 ii=1,n-1 if (ABS(e(ii,ii)).LE.TOL) then DO 3 io=ii+1,n 3 if (ABS(e(io,ii)).GT.TOL) goto 11 IERR=1 write(6,*)ii write(6,*)'tol = ',tol tol=tol*0.50d0 if(tol.gt.1.0d-20)goto 10000 RETURN c 11 CONTINUE DO 4 kk=1,2*n w=e(ii,kk) e(ii,kk)=e(io,kk) 4 e(io,kk)=w ENDIF eii=e(ii,ii) DO 5 jj=ii+1,n e1=e(jj,ii)/eii DO 6 kk=ii+1, 2*n 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e1 5 e(jj,ii)=0.0D0 2 CONTINUE c DO 7 i2=n,2,-1 eii=e(i2,i2) DO 7 j2=i2-1,1,-1 e1=e(j2,i2)/eii DO 9 kk=1, n 9 e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1 7 e(j2,i2)=0.0d0 c DO 10 ii=1,n ei=1.0d0/e(ii,ii) DO 12 jj=1,n 12 ai(ii,jj)=e(ii,jj+n)*ei 10 CONTINUE c RETURN END subroutine vm(so,ns0,ns) implicit none integer*4 ns0,ns,i,j real*8 so(ns0,ns0) do 192 i=1,ns do 192 j=1,ns 192 so(i,j)=0.0d0 return end