program new4dip
      implicit none
      integer*4 NQ,NAT,K,I,L,IERR,j
      real*8 u(3),m(3),a(3,3),ai(3,3),e(3,6),TOL,mo,mn,d,ro,rn,
     1c(3),t(3)
      real*8,allocatable::W(:)
      TOL=1.0d-15
      write(6,*)' F.INP: analyses magnetic dipoles'
      OPEN(17,FILE='F.INP',STATUS='OLD')
      READ(17,*)NQ,NAT,NAT
      allocate(w(NQ))
      DO 1 K=1,NAT
1     READ(17,*)
      write(6,*)'geometry skipped'
      READ(17,*)
      DO 2 L=1,NAT
      DO 2 I=1,NQ
2     READ (17,*) 
      write(6,*)'S-matrix skipped'
      READ(17,*)
      READ(17,*)(w(I),I=1,NQ)
      write(6,*)'Frequencies read'
      do 3 K=1,NQ
      read(17,*,end=99,err=99)u,m
      mo=m(1)*m(1)+m(2)*m(2)+m(3)*m(3)
      ro=u(1)*m(1)+u(2)*m(2)+u(3)*m(3)
      d =u(1)*u(1)+u(2)*u(2)+u(3)*u(3)
      m(1)=ro*u(1)/d
      m(2)=ro*u(2)/d
      m(3)=ro*u(3)/d
      mn=m(1)*m(1)+m(2)*m(2)+m(3)*m(3)
      rn=u(1)*m(1)+u(2)*m(2)+u(3)*m(3)
3     write(6,600)K,w(K),d,mo,mn,ro,rn
600   format(i5,f12.2,5g12.4)

99    write(6,*)'error in dipole reading'
89    close(17)
      end



      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      SUBROUTINE INVs(A,AI,N)
      IMPLICIT none
      integer*4 N,oo,ii,jj,iw,io,kk,i2,j2
      REAL*8 TOL, A(N,N),AI(N,N),w                           
      real*8, allocatable ::E(:,:)
      tol=1.0d-25
      allocate(E(N,2*N))

      DO 1  ii=1,N                                                 
      DO 1  jj=1,N                                                 
      e(ii,jj)=a(ii,jj)
      E(II,JJ+N)=0.0D0                                             
1     if (ii.EQ.jj)e(ii,jj+N)=1.0D0                                

      DO 2 ii=1,N-1                                                
      iw=ii
      if (DABS(e(iw,iw)).LT.TOL) then                                 
       DO 3 io=iw+1,N                                              
       oo=IO
       if (io.GT.N)oo=io-N                                    
3      if (DABS(e(oo,iw)).GE.TOL) goto 11                             
       call report('Inverse cannot be done')
11     CONTINUE                                                       
       DO 4 kk=1, 2*N                                              
       w=e(iw,kk)                                                     
       e(iw,kk)=e(oo,kk)                                              
4      e(oo,kk)=w                                                     
      ENDIF                                                           
      DO 2 jj=ii+1,N                                               
      DO 6 kk=ii+1, 2*N
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii)                    
2     e(jj,ii)=0.0D0                                                  

      DO 7  ii=1, N-1
      i2=N-ii+1                                                    
      DO 7 jj=1,i2-1                                                  
      j2=i2-jj                                                        
      DO 9 kk=1, N                                                 
9     e(j2,kk+N)=e(j2,kk+N)-e(i2,kk+N)*e(j2,i2)/e(i2,i2)
7     e(j2,i2)=0.0D0
 
      DO 10 ii=1,N                                                 
      DO 10 jj=1,N                                                
10    AI(ii,jj)=e(ii,jj+N)/e(ii,ii)                                

      deallocate(E)
      RETURN                                                          
      END