program factor
      implicit none
      integer i,j,k,np,ns,IERR,ii,nl,np1
      real*8,allocatable::x(:),y(:,:),sij(:,:),q(:,:),e(:),ss(:),
     1st(:),sf(:,:)
      real*8 wmin,wmax,xs0,xs1,xs2,norm,a
      character*170,allocatable::sn(:)
      character*80 f,s
      character*10 s10
      wmin=400.0d0
      wmax=3500.0d0
      xs0=0.0d0
      xs1=1.0d0
      xs2=0.0d0
      allocate(sf(2,1))
      sf=0.0d0
c     fast factor analysis
      call getarg(1,f)
      ns=nl(f)
      allocate(sn(ns))
c     load spectral names:
      call lsn(ns,sn,f)
c     get number of points within wmin wmax in first
      np=np1(sn(1),wmin,wmax)
      write(6,*)sn(1)
      if(np.eq.0)stop
      write(6,*)ns,' spectra,',np,' points in first'
      allocate(x(np),y(ns,np),sij(ns,ns),q(ns,ns),e(ns),ss(np),
     1st(np))
      open(9,file=f)
      read(9,80)s
 80   format(a80)
c     read first spectrum
      call ls(x,st,s,wmin,wmax)
c     quick fix baseline
      call qfb(np,x,st)
      a=norm(np,st)
      write(6,599)1,a
 599  format(' Norm of ',i4,':',e12.4)
      st=st/a
      do 21 j=1,np
 21   y(1,j)=st(j)
c     read rest:
      do 2 i=2,ns
      call lsr(ns,x,st,np,i,sn,xs0,xs1,xs2,.false.,1,sf,.false.,1.0d0)
      call qfb(np,x,st)
      a=norm(np,st)
      write(6,599)i,a
      st=st/a
      do 2 j=1,np
 2    y(i,j)=st(j)

      do 3 i=1,ns
      do 3 j=1,ns
      sij(i,j)=0.d0
      do 3 k=1,np
3     sij(i,j)=sij(i,j)+y(i,k)*y(j,k)
      call TRED12(ns,sij,q,e,2,IERR)
      if(IERR.ne.0)then
       write(6,*)'Diagonalization error'
       stop
      endif
      write(6,*)'Eigenvalues:'
      write(6,601)e
601   format(6e12.4)
      do 4 i=1,ns
      write(s10,400)i
 400  format(i10)
      do 5 ii=1,len(s10)
 5    if(s10(ii:ii).ne.' ')goto 6
 6    open(9,file=s10(ii:len(s10))//'.sub.prn')
      ss=0.0d0
      do 8 j=1,np
      do 8 k=1,ns
8     ss(j)=ss(j)+q(k,i)*y(k,j)
      do 7 j=1,np
      if(dabs(ss(j)).lt.1.0d-20)ss(j)=0.0d0
7     write(9,900)x(j),ss(j)
900   format(f12.2,e12.4)
4     close(9)
      end

      function nl(f)
c     number of lines in a file
      implicit none
      character*(*) f
      integer*4 nl,n
      open(9,file=f)
      n=0
1     read(9,*,end=99,err=99)
      n=n+1
      goto 1
99    close(9)
      nl=n
      return
      end

      SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N),Z(N,N),D(*)
      real*8,allocatable::E(:)
      allocate(E(N))
      IF (IEIG .LT. 0) GO TO 110
      DO 100 I = 1, N
C
      DO 100 J = 1, I
      Z(I,J) = A(I,J)
  100 CONTINUE
C
  110 ICODE = IABS(IEIG)
      IF (N .EQ. 1) GO TO 320
C     :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- ::::::::::
      DO 300 II = 2, N
      I = N + 2 - II
      L = I - 1
      H = 0.0D0
      SCALE = 0.0D0
      IF (L .LT. 2) GO TO 130
C     :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) ::::::::::
      DO 120 K = 1, L
  120 SCALE = SCALE + DABS(Z(I,K))
C
      IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140
  130 E(I) = Z(I,L)
      GO TO 290
C
  140 DO 150 K = 1, L
      F = Z(I,K) / SCALE
      Z(I,K) = F
      H = H + F * F
  150 CONTINUE
C
      G = -DSIGN(DSQRT(H),F)
      E(I) = SCALE * G
      H = H - F * G
      Z(I,L) = F - G
      F = 0.0D0
C
      DO 240 J = 1, L
      IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H)
      G = 0.0D0
C     :::::::::: FORM ELEMENT OF A*U ::::::::::
      DO 180 K = 1, J
  180 G = G + Z(J,K) * Z(I,K)
C
      JP1 = J + 1
      IF (L .LT. JP1) GO TO 220
C
      DO 200 K = JP1, L
  200 G = G + Z(K,J) * Z(I,K)
C     :::::::::: FORM ELEMENT OF P ::::::::::
  220 E(J) = G / H
      F = F + E(J) * Z(I,J)
  240 CONTINUE
C
      HH = F / (H + H)
C     :::::::::: FORM REDUCED A ::::::::::
      DO 260 J = 1, L
      F = Z(I,J)
      G = E(J) - HH * F
      E(J) = G
C
      DO 260 K = 1, J
      Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  260 CONTINUE
C
      DO 280 K = 1, L
  280 Z(I,K) = SCALE * Z(I,K)
C
  290 D(I) = H
  300 CONTINUE
C
  320 D(1) = 0.0D0
      E(1) = 0.0D0
C     :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES ::::::::::
      IF (ICODE .NE. 2) GO TO 600
      DO 500 I = 1, N
      L = I - 1
      IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380
C
      DO 360 J = 1, L
      G = 0.0D0
C
      DO 340 K = 1, L
  340 G = G + Z(I,K) * Z(K,J)
C
      DO 360 K = 1, L
      Z(K,J) = Z(K,J) - G * Z(K,I)
  360 CONTINUE
C
  380 D(I) = Z(I,I)
      Z(I,I) = 1.0D0
      IF (L .LT. 1) GO TO 500
C
      DO 400 J = 1, L
      Z(I,J) = 0.0D0
      Z(J,I) = 0.0D0
  400 CONTINUE
C
  500 CONTINUE
C
  620 CALL TQL12 (N,Z,D,IERR,ICODE,E)
      RETURN
C     :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY ::::::::::
  600 DO 610 I=1,N
  610 D(I) = Z(I,I)
      GO TO 620
      END
C
      SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION Z(N,N),D(*),E(*)
      REAL*8 MACHEP
      EPS = 1.0D0
   10 EPS = 0.50D0*EPS
      TOL1 = EPS + 1.0D0
      IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10
      IF(TOL1-EPS.EQ.1.0D0)  EPS = EPS + EPS
      MACHEP = EPS
C
C     MACHEP=16.0D0**(-13)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      F = 0.0D0
      B = 0.0D0
      E(N) = 0.0D0
C
      DO 240 L = 1, N
      J = 0
      H = MACHEP * (DABS(D(L)) + DABS(E(L)))
      IF (B .LT. H) B = H
C     :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT ::::::::::
      DO 110 M = L, N
      IF (DABS(E(M)) .LE. B) GO TO 120
C     :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ::::::::::
  110 CONTINUE
C
  120 IF (M .EQ. L) GO TO 220
  130 IF (J .EQ. 30) GO TO 1000
      J = J + 1
C     :::::::::: FORM SHIFT ::::::::::
      L1 = L + 1
      G = D(L)
      P = (D(L1) - G) / (2.0D0 * E(L))
      R = DSQRT(P*P+1.0D0)
      D(L) = E(L) / (P + DSIGN(R,P))
      H = G - D(L)
C
      DO 140 I = L1, N
  140 D(I) = D(I) - H
C
      F = F + H
C     :::::::::: QL TRANSFORMATION ::::::::::
      P = D(M)
      C = 1.0D0
      S = 0.0D0
      MML = M - L
C     :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
      DO 200 II = 1, MML
      I = M - II
      G = C * E(I)
      H = C * P
      IF (DABS(P) .LT. DABS(E(I))) GO TO 150
      C = E(I) / P
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * P * R
      S = C / R
      C = 1.0D0 / R
      GO TO 160
  150 C = P / E(I)
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * E(I) * R
      S = 1.0D0 / R
      C = C * S
  160 P = C * D(I) - S * G
      D(I+1) = H + S * (C * G + S * D(I))
C     :::::::::: FORM VECTOR ::::::::::
      IF (ICODE .NE. 2) GO TO 200
      DO 180 K = 1, N
      H = Z(K,I+1)
      Z(K,I+1) = S * Z(K,I) + C * H
      Z(K,I) = C * Z(K,I) - S * H
  180 CONTINUE
C
  200 CONTINUE
C
      E(L) = S * P
      D(L) = C * P
      IF (DABS(E(L)) .GT. B) GO TO 130
  220 D(L) = D(L) + F
  240 CONTINUE
C     :::::::::: ORDER EIGENVALUES AND EIGENVECTORS ::::::::::
      DO 300 II = 2, N
      I = II - 1
      K = I
      P = D(I)
C
      DO 260 J = II, N
      IF (D(J) .LE. P) GO TO 260
      K = J
      P = D(J)
  260 CONTINUE
C
      IF (K .EQ. I) GO TO 300
      D(K) = D(I)
      D(I) = P
C
      IF (ICODE .NE. 2) GO TO 300
      DO 280 J = 1, N
      P = Z(J,I)
      Z(J,I) = Z(J,K)
      Z(J,K) = P
  280 CONTINUE
C
  300 CONTINUE
C
      GO TO 1001
C     :::::::::: SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ::::::::::
 1000 IERR = L
 1001 RETURN
      END

      subroutine lsn(ns,sn,s)
      implicit none
      integer*4 j,ns
      character*(*) s
      character*170 sn(ns)
      open(9,file=s)
      do 1 j=1,ns
1     read(9,9000,end=99,err=99)sn(j)
9000  format(a170)
99    close(9)
      write(6,*)'names loaded'
      return
      end

      subroutine lsr(ns,xs,ys,np,js,sn,xs0,xs1,xs2,lstable,nt,sf,
     1lcon,g)
c     loads in the full spectrum, adjust number of points to np
c     with the scaled xs
      implicit none
      integer*4 np,i,js,npn,j,nt,npall,ns
      real*8 xs(*),ys(*),xi,xj,xjp,xs0,xs1,xs2,sf(2,nt),xjf,g,
     1ef,pi
      real*8 ,allocatable::sx(:),sy(:),sp(:)
      character*170 sn(ns)
      logical lstable,lcon
      pi=4.0d0*datan(1.0d0)
      npn=npall(sn(js))
      allocate(sx(npn),sy(npn),sp(npn))
c     loads in the full spectrum:
      call lsf(npn,sx,sy,sn(js))
c     convolute
      if(lcon)then
       do 3 i=1,npn
       sp(i)=0.0d0
       do 31 j=1,npn-1
       ef=((sx(i)-sx(j))/g)**2
31     if(ef.lt.20.0d0)sp(i)=sp(i)+exp(-ef)*sy(j)*dabs(sx(j)-sx(j+1))
3      sp(i)=sp(i)/dsqrt(pi)/g
       sy=sp
      endif
      do 1 i=1,np
      xi=xs(i)
      ys(i)=0.0d0
      do 2 j=1,npn-1
c     scaled values:
      if(lstable)then
       xj =xjf(sx(j  ),nt,sf)
       xjp=xjf(sx(j+1),nt,sf)
      else
       xj =xs0+xs1*sx(j  )+xs2*sx(j  )*2
       xjp=xs0+xs1*sx(j+1)+xs2*sx(j+1)*2
      endif
      if(xj.le.xi.and.xjp.ge.xi)then
       ys(i)=sy(j)+(sy(j+1)-sy(j))*(xi-xj)/(xjp-xj)
       goto 1
      endif
 2    continue
 1    continue
      return
      end

      function xjf(xj,nt,sf)
      implicit none
      integer*4 nt,i
      real*8 xjf,xj,sf(2,nt),xe1,xe2,xc1,xc2
      do 1 i=1,nt-1
      xe1=sf(1,i )
      xe2=sf(1,i+1)
      xc1=sf(2,i )
      xc2=sf(2,i+1)
c     xe1,xe2 ... experimental values, xc1,xc2 ... corresponding calc.
      if((xj.ge.xc1.and.xj.le.xc2).or.
     1   (xj.le.xc1.and.xj.ge.xc2))then
       xjf=xe1+(xj-xc1)*(xe2-xe1)/(xc2-xc1)
       return
      endif
 1    continue
      if(dabs(xj-sf(2,nt)).lt.dabs(xj-sf(2,1)))then
       xjf=sf(1,nt)+xj-sf(2,nt)
      else
       xjf=sf(1,1)+xj-sf(2,1)
      endif
      return
      end

      function npall(s)
      implicit none
      integer*4 npall,n
      character*(*) s
      n=0
      open(9,file=s)
101   read(9,*,end=99,err=99)
      n=n+1
      goto 101
99    close(9)
      npall=n
      return
      end

      function np1(s,wmin,wmax)
      implicit none
      integer*4 np1,n
      real*8 wmin,wmax,x
      character*(*) s
      open(9,file=s)
      n=0
101   read(9,*,end=99,err=99)x
      if(x.ge.wmin.and.x.le.wmax)n=n+1
      goto 101
99    close(9)
      np1=n
      return
      end

      subroutine ls(xs,ys,name,wmin,wmax)
c     loads in the spectrum  within wmin wmax    
      implicit none
      integer*4 i
      real*8 xs(*),ys(*),wmin,wmax,x,y
      character*(*) name
      i=0
      open(9,file=name)
 11   read(9,*,end=99,err=99)x,y
      if(wmin.le.x.and.x.le.wmax)then
       i=i+1       
       xs(i)=x
       ys(i)=y
      endif
      goto 11
  99  close(9)
      return
      end

      function norm(np,s)
      implicit none
      integer*4 i,np
      real*8 norm,s(*)
      norm=0.0d0
      do 1 i=1,np
   1  norm=norm+s(i)**2
      norm=dsqrt(norm)
      return
      end

      subroutine lsf(n,xs,ys,name)
c     loads in the full spectrum      
      implicit none
      integer*4 i,n
      real*8 xs(*),ys(*)
      character*(*) name
      open(9,file=name)
      do 11 i=1,n
 11   read(9,*)xs(i),ys(i)
      close(9)
      return
      end

      subroutine qfb(np,x,y)
      implicit none
      integer*4 np,i
      real*8 x(np),y(np),x1,x2,y1,y2,xa,xb,aa,bb
      xa=1900.0d0
      xb=2600.0d0
      aa=dabs(x(1)-xa)
      bb=dabs(x(1)-xb)
      x1=x(1)
      x2=x(1)
      y1=y(1)
      y2=y(2)
      do 1 i=1,np
      if(dabs(x(i)-xa).lt.aa)then
       aa=dabs(x(i)-xa)
       x1=x(i)
       y1=y(i)
      endif
      if(dabs(x(i)-xb).lt.bb)then
       bb=dabs(x(i)-xb)
       x2=x(i)
       y2=y(i)
      endif
1     continue
      do 2 i=1,np
2     y(i)=y(i)-y1-(y2-y1)/(x2-x1)*(x(i)-x1)
      return
      end