      PROGRAM swath
      IMPLICIT NONE
c----------------------------------------------
c DICHIARAZIONI E DIMENSIONAMENTI
c----------------------------------------------

      INTEGER i,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,j
      INTEGER nbuff,ntot,INFO,k
      INTEGER neq,ntime,iradia_diffra,nperiodi
      INTEGER IPVT(3000)

      REAL*8 l,d,t,w,a,r,teta,alfa,dx,dy,dpsi,pi,rho,grav,alamb
      REAL*8 x1(3000),x2(3000),y1(3000),y2(3000),xmed(3000),ymed(3000),
     $       dist(3000),nx(3000),ny(3000)
      REAL*8 C(3000,3000,2),AA(3000,3000)
      real*8 yin(2*3000),yout(2*3000),dyindt(2*3000)
      real*8 fi(3000),fiold(3000),dfidn(3000),press(3000),dfidt(3000)
      real*8 time,t_sub,omega,amp_sway,amp_heave,
     $       amp_wavemaker,amu,amp_wave,eta_indist
      real*8 fx,fy,fxsx,fysx,fxdx,fydx,fxond,fyond,xenc,yenc,dt,periodo
      real*8 xx(500),yy(500)

c----------------------------------------------
c LETTURA DATI DELLA VASCA
c----------------------------------------------
c
      OPEN(unit=20,file='Dati_mono_hull.dat',
     $access='sequential',form='formatted')
	  Read(20,*) l   !LUNGHEZZA DEL DOMINIO
	  Read(20,*) d   !PROFONDITA
	  Read(20,*) n1,n2,n3,n4,n5   
      Read(20,*)amp_sway,amp_heave,amp_wavemaker,amp_wave,alamb
      Read(20,*)iradia_diffra    !0=RADIAZIONE    1= DIFFRAZIONE
	  Read(20,*)nperiodi
	  Read(20,*)amu   
      close(20)

c----------------------------------------------
c LETTURA DATI DELLA SEZIONE TRASVERSALE
c----------------------------------------------
c
      OPEN(unit=22,file='Rilevato.dat',
     $access='sequential',form='formatted')
      do i=1,n3+1
          read(22,*)xx(i),yy(i)
      end do
      close(22)

c----------------------------------------------
c APERTURA FILES DI OUTPUT
c----------------------------------------------
c
      open(15,file='forze.dat',status='unknown',
     >     form='formatted',access='sequential')
      open(88,file='elev_wavem.dat',status='unknown',
     >     form='formatted',access='sequential')

      open(2,file='anima.plt',status='unknown',   !file di animazione grafica su Tecplot
     >     form='formatted',access='sequential')
      write(2,301)
      write(2,302)
  301 format('title="animazione"') 
  302 format('variables="x","y","px","py"')  
   
c----------------------------------------------
c DATI GENERALI
c----------------------------------------------
c
      pi=acos(-1.)
      grav=9.8065d0
      rho=1000.d0
      amu=amu*pi/180.d0

      omega=sqrt(2.*pi*grav*tanh(2*pi*d/alamb)/alamb)
      write(*,*)omega
      pause
      periodo=2*pi/omega
      dt=periodo/40.
      ntime=40*nperiodi
    
c----------------------------------------------
c GRIGLIA DI CALCOLO
c----------------------------------------------
c
c PARETE SINISTRA
c      
      dy=abs(d)/n1

      DO i=1,n1
        x1(i)=0
        y1(i)=-d+(i-1)*dy
        x2(i)=0
        y2(i)=-d+i*dy
      END DO

c SUPERFICIE LIBERA A SINISTRA DELLA CARENA      
c
      dx=(l/2.d0+xx(1))/n2
      nbuff=n1

      DO i=nbuff+1,nbuff+n2
	    x1(i)=(i-nbuff-1)*dx
	    y1(i)=0
	    x2(i)=(i-nbuff)*dx
	    y2(i)=0
      END DO

c IMESH DI CARENA
c
      DO i=1,n3
        x1(i+n1+n2)= xx(i)+l/2
		y1(i+n1+n2)= yy(i)
        x2(i+n1+n2)= xx(i+1)+l/2
		y2(i+n1+n2)= yy(i+1)
      END DO

c SUPERFICIE LIBERA A DESTRA DELLA CARENA
c     
      dx=(l/2-xx(n3+1))/n4
      nbuff=n1+n2+n3

	  DO i=nbuff+1,nbuff+n4
	    x1(i)=l/2+xx(n3+1)+(i-nbuff-1)*dx
	    y1(i)=0
	    x2(i)=l/2+xx(n3+1)+(i-nbuff)*dx
	    y2(i)=0
      END DO

c PARETE DESTRA
c      
      dy=abs(d)/n5
      nbuff=n1+n2+n3+n4

	  DO i=nbuff+1,nbuff+n5
	    x1(i)=l
	    y1(i)=-(i-nbuff-1)*dy
	    x2(i)=l
	    y2(i)=-(i-nbuff)*dy
      END DO

c NUMERO TOTALE DI ELEMENTI
c
      ntot=nbuff+n5

c PUNTI MEDI, SEGMENTI E VERSORI
c
      DO i=1,ntot
        xmed(i)= (x1(i)+x2(i))/2
        ymed(i)= (y1(i)+y2(i))/2
	    dist(i)= sqrt((x2(i)-x1(i))**2+(y2(i)-y1(i))**2)
	    nx(i) =  (y2(i)-y1(i))/dist(i)
        ny(i) = -(x2(i)-x1(i))/dist(i)     
      END DO

c----------------------------------------------
c CALCOLO DEI KERNELS - SECONDA IDENTITA' DI GREEN - METODO DIRETTO
c----------------------------------------------
c
      call integrali(x1,xmed,x2,y1,ymed,y2,nx,ny,dist,
     $               C,d,pi,ntot)

c----------------------------------------------
c SET-UP DELLA MATRICE 
c----------------------------------------------
c
	call matrix(C,AA,PI,ntot,n1,n2,n3,n4,n5)

c----------------------------------------------
c INVERSIONE DEFINITIVA DELLA MATRICE
c----------------------------------------------
c
	call DGEFA (AA,3000,ntot,IPVT,INFO)

c----------------------------------------------
C NUMERO DI EQUAZIONI DA INTEGRARE NEL TEMPO
c----------------------------------------------
c
	neq=2*(n2+n4)

c----------------------------------------------
C INIZIALIZZAZIONE DELLE VARIABILI
c----------------------------------------------
c
      time=0.d0

      do i=1,neq
        yin(i)=0.d0
      end do

      do i=1,ntot
        fi(i)=0.d0
        fiold(i)=0.d0
        dfidn(i)=0.d0
      enddo

c----------------------------------------------
c Stampa per animazione TECPLOT
c----------------------------------------------
      write(2,305)int(time),ntot     
      do i=1,ntot
        write(2,1002)xmed(i)+amp_sway *cos(omega*time),
     $               ymed(i)+amp_heave*cos(omega*time),0.,0.
      end do

c----------------------------------------------
c LOOP SUL TEMPO
c----------------------------------------------
c 
      do k=1,ntime+1
        write(*,*)
        write(*,*)'******   ntime=',k,'   t=',time/periodo,'   ******'

c------------------
c Calcolo delle derivate dyindt delle equazioni da avanzare nel tempo
c------------------
        call derivs(time,yin,dyindt,AA,C,IPVT,
     $              ntot,n1,n2,n3,n4,n5,
     $              omega,amp_sway,amp_heave,amp_wavemaker,amp_wave,
     $              pi,grav,fi,
     $              time,alamb,iradia_diffra,xmed,ymed,nx,ny,d,amu)

c------------------
c Inizializzo le forze sulla carena complessiva
c------------------
        fx=0.d0
	    fy=0.d0
c------------------
c Inizializzo le forze sull'ondogeno
c------------------
	    fxond=0.d0
	    fyond=0.d0

c------------------
c Scelta del caso RADIA O DIFFRA per il calcolo delle pressioni e forze sulla carena 
c------------------
        if (iradia_diffra.eq.0) then   !radiazione o diffrazione con ondogeno
          do i = n1+n2+1,n1+n2+n3                                        
            dfidt(i)=(fi(i)-fiold(i))/DT
	        press(i)=-rho*dfidt(i)
            fx=fx-press(i)*dist(i)*nx(i)
	        fy=fy-press(i)*dist(i)*ny(i)
            fiold(i)=fi(i)
          enddo
        else   !diffrazione col metodo della scomposizione del potenziale       
          do i = n1+n2+1,n1+n2+n3         
            dfidt(i)=(fi(i)-fiold(i))/DT
            xenc= xmed(i)*sin(amu)
            yenc= ymed(i)
	        press(i)=-rho*dfidt(i)+
     $                 rho*amp_wave*grav*
     $                (cosh(2*pi/alamb*(d+yenc))/
     $                 cosh(2*pi/alamb*d))*
     $                 cos(2*pi/alamb*xenc-omega*time)
            fx=fx-press(i)*dist(i)*nx(i)
	        fy=fy-press(i)*dist(i)*ny(i)
            fiold(i)=fi(i)
	      enddo
	    endif	

c------------------
c Pessioni e forze sull' ondogeno
c------------------
         do i = 1,n1
           dfidt(i)=(fi(i)-fiold(i))/DT
	       press(i)=-rho*dfidt(i)
	       fxond=fxond-press(i)*dist(i)*nx(i)
	       fyond=fyond-press(i)*dist(i)*ny(i)
	       fiold(i)=fi(i)
	    enddo

        write(15,1100)time/periodo,fx,fy

c------------------
c Schema di RK4 per l'avanzamento nel tempo
c------------------
        call rk4(yin,dyindt,neq,time,dt,yout,AA,C,IPVT,
     $           ntot,n1,n2,n3,n4,n5,
     $           omega,amp_sway,amp_heave,amp_wavemaker,amp_wave,
     $           pi,grav,fi,
     $           alamb,iradia_diffra,xmed,ymed,nx,ny,d,amu)

c------------------
c Aggiornamento delle variabili prima di ripartire con il nuovo time step
c------------------
        do i = 1,neq
          yin(i)=yout(i)
        end do
      
c------------------
c File di grafica per Tecplot
c------------------
        write(2,305)k,ntot

        if (iradia_diffra.eq.0) then   !radiazione o diffrazione con ondogeno
          do i=1,n1
            write(2,1002)xmed(i),ymed(i),0.,0.
          end do
          do i=1,n2
            write(2,1002)xmed(n1+i),yout(i),0.,0.
          end do
          do i= n1+n2+1,n1+n2+n3
             write(2,1002)xmed(i)+amp_sway*cos(omega*(time+dt)),
     $                    ymed(i)+amp_heave*cos(omega*(time+dt)),
     $                    -press(i)*nx(i),-press(i)*ny(i)
          end do
          do i=1,n4
            write(2,1002)xmed(n1+n2+n3+i),yout(i+n2),0.,0.
          end do
          do i=1,n5
             write(2,1002)xmed(n1+n2+n3+n4+i),
     $                    ymed(n1+n2+n3+n4+i),0.,0.
          end do
          write(88,*)time/periodo,yout(1),
     $               amp_wavemaker*cos(omega*(time+dt))
        else        !diffrazione col metodo della scomposizione del potenziale
          do i=1,n1
             write(2,1002)xmed(i),ymed(i),0.,0.
          end do
          do i=1,n2
            xenc = xmed(n1+i)*sin(amu)
	        eta_indist=amp_wave*cos(2*pi/alamb*xenc-omega*(time+dt))
            write(2,1002)xmed(n1+i),yout(i)+eta_indist,0.,0.
c               write(2,1002)xmed(n1+i),yout(i),0.,0.
          end do
          do i= n1+n2+1,n1+n2+n3
             write(2,1002)xmed(i)+amp_sway*cos(omega*(time+dt)),
     $                    ymed(i)+amp_heave*cos(omega*(time+dt)),
     $                    -press(i)*nx(i),-press(i)*ny(i)
          end do
          do i=1,n4
            xenc = xmed(n1+n2+n3+i)*sin(amu)
	        eta_indist=amp_wave*cos(2*pi/alamb*xenc-omega*(time+dt))
            write(2,1002)xmed(n1+n2+n3+i),
     $                      yout(i+n2)+eta_indist,0.,0.
c     $                      yout(i+n2),0.,0.
          end do
          do i=1,n5
             write(2,1002)xmed(n1+n2+n3+n4+i),
     $                    ymed(n1+n2+n3+n4+i),0.,0.
          end do
        end if

c------------------
c Aggiornamento del tempo
c------------------
        time=time+dt
      end do
c--------------------------------------------------------------
      close (15)
      close (88)

c-------------------------------------------
  305 format('zone t="',i6,1x,'",','i=',i7,1x,'f=point')
 1002 format (6(f15.5,1x))
 1100 format (7(f15.5,1x))
      stop
	end

c***********************************************************************
      subroutine integrali(x1,x,x2,y1,y,y2,anx,any,dim,
     >                  c,d,pigr,ntot)

      implicit real*8(a-h,o-z)
      parameter (nmax=3000)

      dimension x1(NMAX),x(NMAX),x2(NMAX),
     >          y1(NMAX),y(NMAX),y2(NMAX),
     >          anx(NMAX),any(NMAX),dim(NMAX)
      dimension c(nmax,nmax,2)
c
      do i=1,ntot
         do j=1,ntot
            eps=1.0e-12
c
c alla Jaswon-Symm (integrali in senso antiorario normale verso l'interno)
            dist1=sqrt((x2(j)-x(i))**2+(y2(j)-y(i))**2)
            dist2=sqrt((x1(j)-x(i))**2+(y1(j)-y(i))**2)
            arg_psi=(dist1**2+dist2**2-dim(j)**2)/2.d0/dist1/dist2
            arg_theta=(dist1**2+dim(j)**2-dist2**2)/2./dist1/dim(j)
            if(abs(arg_psi).gt.1.d0) arg_psi=sign(1.d0,arg_psi)
            if(abs(arg_theta).gt.1.d0) arg_theta=sign(1.d0,arg_theta)
            psi=acos(arg_psi)
            theta=acos(arg_theta)
            scalar=anx(j)*(x(i)-x(j))+any(j)*(y(i)-y(j))
        
            if(i.eq.j) then
               c(i,j,1)=-pigr
            elseif(scalar.eq.0.d0) then
               c(i,j,1)=0.
            elseif(scalar.gt.0.d0) then
                c(i,j,1)=psi
            elseif(scalar.lt.0.d0) then
               c(i,j,1)=-psi
            end if
         
            if(i.eq.j) then
               c(i,j,2)=dim(j)*(log(dim(j)/2.d0)-1.d0)
            elseif(theta.le.eps) then
               c(i,j,2)=dist1*(log(dist1)-1.d0)-dist2*(log(dist2)-1.d0)
            elseif(pigr-theta.le.eps) then
               c(i,j,2)=-dist1*(log(dist1)-1.d0)+dist2*(log(dist2)-1.d0)
            else
               c(i,j,2)=dist1*cos(theta)*(log(dist1)-log(dist2))+
     >             dim(j)*(log(dist2)-1.d0)+dist1*psi*sin(theta)
            end if
        
c pannello simmetrico rispetto a y=-d

            dist1=sqrt((x1(j)-x(i))**2+(-2.*d-y1(j)-y(i))**2)  
            dist2=sqrt((x2(j)-x(i))**2+(-2.*d-y2(j)-y(i))**2)
            arg_psi=(dist1**2+dist2**2-dim(j)**2)/2./dist1/dist2
            arg_theta=(dist1**2+dim(j)**2-dist2**2)/2.d0/dist1/dim(j)
            if(abs(arg_psi).gt.1.d0) arg_psi=sign(1.d0,arg_psi)
            if(abs(arg_theta).gt.1.d0) arg_theta=sign(1.d0,arg_theta)
            psi=acos(arg_psi)
            theta=acos(arg_theta)
            scalar=anx(j)*(x(i)-x(j))-any(j)*(y(i)+2.*d+y(j))
 
            if(scalar.eq.0.d0) then
               c(i,j,1)=c(i,j,1)
            elseif(scalar.gt.0.d0) then
               c(i,j,1)=c(i,j,1)+psi
            elseif(scalar.lt.0.d0) then
               c(i,j,1)=c(i,j,1)-psi
            end if
 
            if(theta.le.eps) then
               c(i,j,2)=c(i,j,2)+dist1*(log(dist1)-1.d0)-
     >		 dist2*(log(dist2)-1.d0)
            elseif(pigr-theta.le.eps) then
               c(i,j,2)=c(i,j,2)-dist1*(log(dist1)-1.d0)+
     >		 dist2*(log(dist2)-1.d0)
            else
               c(i,j,2)=c(i,j,2)+dist1*cos(theta)*(log(dist1)-
     >		 log(dist2))+
     >              dim(j)*(log(dist2)-1.d0)+dist1*psi*sin(theta)
            end if
      
          end do
      end do
c
      return
      end 
c***********************************************************************
      subroutine matrix(C,A,
     $PIGR,ntot,n1,n2,n3,n4,n5)
	      implicit real*8(a-h,o-z)
      parameter (nmax=3000)

      dimension a(NMAX,nmax)
      dimension c(nmax,nmax,2)

       do i=1,ntot
 
          do j=1,n1   
             a(i,j)=c(i,j,1)    !Neumann
             end do

          do j=n1+1,n1+n2
             a(i,j)=c(i,j,2)     !Dirichlet
		end do
c       
          do j=n1+n2+1,n1+n2+n3
             a(i,j)=c(i,j,1)    !Neumann
          end do
c
          do j=n1+n2+n3+1,n1+n2+n3+n4
             a(i,j)=c(i,j,2)     !Dirichlet
          end do
c
          do j=n1+n2+n3+n4+1,n1+n2+n3+n4+n5
             a(i,j)=c(i,j,1)    !Neumann
          end do
c
        end do
      return
	end
c***********************************************************************
      subroutine termini_noti(C,fi,DFIDN,B,
     $PIGR,ntot,n1,n2,n3,n4,n5)
	      implicit real*8(a-h,o-z)
      parameter (nmax=3000)

      dimension B(NMAX),fi(nmax),dfidn(nmax)
      dimension c(nmax,nmax,2)


	do i=1,ntot
          b(i)=0.d0
          do j=1,n1
             b(i)=b(i)-c(i,j,2)*dfidn(j)
          end do
c              
          do j=n1+1,n1+n2
             b(i)=b(i)-c(i,j,1)*fi(j)
          end do
c       
          do j=n1+n2+1,n1+n2+n3
             b(i)=b(i)-c(i,j,2)*dfidn(j)
          end do
c
          do j=n1+n2+n3+1,n1+n2+n3+n4
             b(i)=b(i)-c(i,j,1)*fi(j)
          end do
c
          do j=n1+n2+n3+n4+1,n1+n2+n3+n4+n5
             b(i)=b(i)-c(i,j,2)*dfidn(j)
          end do
c
      end do
c  
       return
       end 
c****************************************************************************

c***********************************************************************
      SUBROUTINE DGEFA (A,LDA,N,IPVT,INFO)
      INTEGER LDA,N,IPVT(1),INFO
      DOUBLE PRECISION A(LDA,1)
C
C     DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
C
C     DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
C     (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE MATRIX TO BE FACTORED.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C     ON RETURN
C
C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C                WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(N)
C                AN INTEGER VECTOR OF PIVOT INDICES.
C
C        INFO    INTEGER
C                = 0  NORMAL VALUE.
C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C                     INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
C                     IF CALLED.  USE  RCOND  IN DGECO FOR A RELIABLE
C                     INDICATION OF SINGULARITY.
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DSCAL,IDAMAX
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION T
      INTEGER IDAMAX,J,K,KP1,L,NM1
C
C
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      INFO = 0
      NM1 = N - 1
      IF (NM1 .LT. 1) GO TO 70
      DO 60 K = 1, NM1
         KP1 = K + 1
C
C        FIND L = PIVOT INDEX
C
         L = IDAMAX(N-K+1,A(K,K),1) + K - 1
         IPVT(K) = L
C
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
         IF (A(L,K) .EQ. 0.0D0) GO TO 40
C
C           INTERCHANGE IF NECESSARY
C
            IF (L .EQ. K) GO TO 10
               T = A(L,K)
               A(L,K) = A(K,K)
               A(K,K) = T
   10       CONTINUE
C
C           COMPUTE MULTIPLIERS
C
            T = -1.0D0/A(K,K)
            CALL DSCAL(N-K,T,A(K+1,K),1)
C
C           ROW ELIMINATION WITH COLUMN INDEXING
C
            DO 30 J = KP1, N
               T = A(L,J)
               IF (L .EQ. K) GO TO 20
                  A(L,J) = A(K,J)
                  A(K,J) = T
   20          CONTINUE
               CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
   30       CONTINUE
         GO TO 50
   40    CONTINUE
            INFO = K
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      IPVT(N) = N
      IF (A(N,N) .EQ. 0.0D0) INFO = N
      RETURN
      END
      SUBROUTINE DGESL (A,LDA,N,IPVT,B,JOB)
      INTEGER LDA,N,IPVT(1),JOB
      DOUBLE PRECISION A(LDA,1),B(1)
C
C     DGESL SOLVES THE DOUBLE PRECISION SYSTEM
C     A * X = B  OR  TRANS(A) * X = B
C     USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE OUTPUT FROM DGECO OR DGEFA.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C        IPVT    INTEGER(N)
C                THE PIVOT VECTOR FROM DGECO OR DGEFA.
C
C        B       DOUBLE PRECISION(N)
C                THE RIGHT HAND SIDE VECTOR.
C
C        JOB     INTEGER
C                = 0         TO SOLVE  A*X = B ,
C                = NONZERO   TO SOLVE  TRANS(A)*X = B  WHERE
C                            TRANS(A)  IS THE TRANSPOSE.
C
C     ON RETURN
C
C        B       THE SOLUTION VECTOR  X .
C
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
C        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE
C        CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0
C        OR DGEFA HAS SET INFO .EQ. 0 .
C
C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
C     WITH  P  COLUMNS
C           CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
C           IF (RCOND IS TOO SMALL) GO TO ...
C           DO 10 J = 1, P
C              CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
C        10 CONTINUE
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DDOT
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION DDOT,T
      INTEGER K,KB,L,NM1
C
C
      NM1 = N - 1
      IF (JOB .NE. 0) GO TO 50
C
C        JOB = 0 , SOLVE  A * X = B
C        FIRST SOLVE  L*Y = B
C
         IF (NM1 .LT. 1) GO TO 30
         DO 20 K = 1, NM1
            L = IPVT(K)
            T = B(L)
            IF (L .EQ. K) GO TO 10
               B(L) = B(K)
               B(K) = T
   10       CONTINUE
            CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
   20    CONTINUE
   30    CONTINUE
C
C        NOW SOLVE  U*X = Y
C
         DO 40 KB = 1, N
            K = N + 1 - KB
            B(K) = B(K)/A(K,K)
            T = -B(K)
            CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
   40    CONTINUE
      GO TO 100
   50 CONTINUE
C
C        JOB = NONZERO, SOLVE  TRANS(A) * X = B
C        FIRST SOLVE  TRANS(U)*Y = B
C
         DO 60 K = 1, N
            T = DDOT(K-1,A(1,K),1,B(1),1)
            B(K) = (B(K) - T)/A(K,K)
   60    CONTINUE
C
C        NOW SOLVE TRANS(L)*X = Y
C
         IF (NM1 .LT. 1) GO TO 90
         DO 80 KB = 1, NM1
            K = N - KB
            B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
            L = IPVT(K)
            IF (L .EQ. K) GO TO 70
               T = B(L)
               B(L) = B(K)
               B(K) = T
   70       CONTINUE
   80    CONTINUE
   90    CONTINUE
  100 CONTINUE
      RETURN
      END
      INTEGER FUNCTION IDAMAX(N,DX,INCX)
C
C     FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DMAX
      INTEGER I,INCX,IX,N
C
      IDAMAX = 0
      IF( N .LT. 1 ) RETURN
      IDAMAX = 1
      IF(N.EQ.1)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      IX = 1
      DMAX = DABS(DX(1))
      IX = IX + INCX
      DO 10 I = 2,N
         IF(DABS(DX(IX)).LE.DMAX) GO TO 5
         IDAMAX = I
         DMAX = DABS(DX(IX))
    5    IX = IX + INCX
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
   20 DMAX = DABS(DX(1))
      DO 30 I = 2,N
         IF(DABS(DX(I)).LE.DMAX) GO TO 30
         IDAMAX = I
         DMAX = DABS(DX(I))
   30 CONTINUE
      RETURN
      END
      SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
C
C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DA
      INTEGER I,INCX,INCY,IXIY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF (DA .EQ. 0.0D0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DY(IY) + DA*DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,4)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DY(I) + DA*DX(I)
   30 CONTINUE
      IF( N .LT. 4 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
        DY(I) = DY(I) + DA*DX(I)
        DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
        DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
        DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
   50 CONTINUE
      RETURN
      END      
      
      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      DDOT = 0.0D0
      DTEMP = 0.0D0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DTEMP + DX(IX)*DY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      DDOT = DTEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DTEMP + DX(I)*DY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
     *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
   50 CONTINUE
   60 DDOT = DTEMP
      RETURN
      END
      SUBROUTINE  DSCAL(N,DA,DX,INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DA,DX(1)
      INTEGER I,INCX,M,MP1,N,NINCX
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
        DX(I) = DA*DX(I)
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DX(I) = DA*DX(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DX(I) = DA*DX(I)
        DX(I + 1) = DA*DX(I + 1)
        DX(I + 2) = DA*DX(I + 2)
        DX(I + 3) = DA*DX(I + 3)
        DX(I + 4) = DA*DX(I + 4)
   50 CONTINUE
      RETURN
      END
c***********************************************************************
	
c***********************************************************************
      subroutine derivs(t,yin,dyindt,AA,C,IPVT,ntot,n1,n2,n3,n4,
     $                  n5,
     $                  omega,amp_sway,amp_heave,amp_wavemaker,amp_wave,
     $                  pigr,grav,fi,t_sub,
     $                  alamb,iradia_diffra,xmed,ymed,anx,any,d,amu)

      implicit real*8(a-h,o-z)
      parameter (nmax=3000,nmax_eq=2*nmax)
      dimension eta(nmax),fi(nmax),dfidn(nmax),
     $          xmed(nmax),ymed(nmax),anx(nmax),any(nmax)
      dimension AA(nmax,nmax),C(nmax,nmax,2),IPVT(nmax),b(nmax)
      dimension yin(NMAX_EQ),dyindt(NMAX_EQ),u(nmax),v(nmax)
c


c----------------------------------------------
c AGGIORNAMENTO DEL POTENZIALE E DELL'ELEVAZIONE DELLA SUPERFICIE LIBERA DI SINISTRA
c----------------------------------------------
      do i=1,n2
        j=i
        eta(j)=yin(j)     
        fi(n1+i)=yin(i+n2+n4)
      end do

c----------------------------------------------
c AGGIORNAMENTO DEL POTENZIALE E DELL'ELEVAZIONE DELLA SUPERFICIE LIBERA DI DESTRA
c----------------------------------------------
      do i=1,n4
        j=i+n2
        eta(j)=yin(j)     
        fi(n1+n2+n3+i)=yin(j+n2+n4)
      end do

c----------------------------------------------
c AGGIORNAMENTO DELLA CONDIZIONE AL CONTORNO SULLA PARETE ONDOGENO
c----------------------------------------------
      if (iradia_diffra.eq.0) then   !radiazione o diffrazione con ondogeno
        do i=1,n1                              
          dfidn(i)=-amp_wavemaker*omega*sin(omega*t)
	    enddo
      else                !diffrazione col metodo della scomposizione del potenziale
        do i=1,n1                                     
          dfidn(i)=0.
        enddo
      end if

c----------------------------------------------
c AGGIORNAMENTO DELLA CONDIZIONE AL CONTORNO SULLA PARETE DI DESTRA
c----------------------------------------------
      do i=1+n1+n2+n3+n4,ntot                                     
	    dfidn(i)=0.d0
      enddo

c----------------------------------------------
c AGGIORNAMENTO DELLA CONDIZIONE AL CONTORNO SULLA CARENA
c----------------------------------------------
      do i=n1+n2+1 , n1+n2+n3                                     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if (iradia_diffra.eq.0) then     !radiazione o diffrazione con ondogeno
	      dfidn(i)=-amp_heave*omega*sin(omega*t_sub)*any(i) !sussulto
     $             -amp_sway *omega*sin(omega*t_sub)*anx(i) !squasso	
        else                             !diffrazione col metodo della scomposizione del potenziale
          xenc= xmed(i)*sin(amu)
          yenc= ymed(i)
	      u(i)=.5*amp_wave*omega*cosh(2*pigr/alamb*(d+yenc))/
     $	       sinh(2*pigr/alamb*d)*cos(2*pigr/alamb*xenc-omega*t_sub)
	      v(i)=.5*amp_wave*omega*sinh(2*pigr/alamb*(d+yenc))/
     $	       sinh(2*pigr/alamb*d)*sin(2*pigr/alamb*xenc-omega*t_sub)
          uenc= u(i)*sin(amu)
          venc= v(i)
	      dfidn(i)=-(uenc*anx(i)+venc*any(i))
	    endif
      enddo

c----------------------------------------------
c CALCOLO DEI TERMINI NOTI DEL SISTEMA LINEARE ASSOCIATO ALLA SECONDA IDENTITA' Di GREEN - METODO DIRETTO
c----------------------------------------------
      call termini_noti(C,fi,DFIDN,B,
     $                  PIGR,ntot,n1,n2,n3,n4,n5)

c----------------------------------------------
c CALCOLO DELLA SOLUZIONE DEL SISTEMA LINEARE ASSOCIATO ALLA SECONDA IDENTITA' Di GREEN - METODO DIRETTO
c----------------------------------------------	
      call dgesl(AA,NMAX,ntot,IPVT,B,0)

c----------------------------------------------
c IDENTIFICAZIONE DEL VETTORE SOLUZIONE NEI CONFRONTI DELLE CONDIZIONI AL CONTORNO
c----------------------------------------------	
      do i=1,n1                                      
        fi(i)=b(i)
      enddo
      do i=n1+1,n1+n2                                     
        dfidn(i)=b(i)
      enddo
      do i=n1+n2+1,n1+n2+n3                                 
        fi(i)=b(i)
      enddo
      do i=n1+n2+n3+1,n1+n2+n3+n4                               
        dfidn(i)=b(i)
      enddo
      do i=n1+n2+n3+n4+1,n1+n2+n3+n4+n5                             
        fi(i)=b(i)
      enddo

c----------------------------------------------
c CALCOLO DELLE DERIVATE DELLA SUPERFICIE LIBERA DI SINISTRA PER L'AVANZAMENTO NEL TEMPO
c----------------------------------------------	
      do i=1,n2
        j=i+n1                                      
        dyindt(i)     = -dfidn(j)          
        dyindt(i+n2+n4)=-grav*eta(i)        
      end do
c----------------------------------------------
c CALCOLO DELLE DERIVATE DELLA SUPERFICIE LIBERA DI DESTRA PER L'AVANZAMENTO NEL TEMPO
c----------------------------------------------	
      do i=1,n4                                      
        j=i+n1+n2+n3                                      
        dyindt(i+n2)     = -dfidn(j)        
        dyindt(i+2*n2+n4)=-grav*eta(i+n2)        
      end do
      
c
      return
      end
c*************************************************************************
c
c*************************************************************************
      SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,AA,C,IPVT,
     $               ntot,n1,n2,n3,n4,n5,
     $               omega,amp_sway,amp_heave,amp_wavemaker,amp_wave,
     $               pigr,grav,fi,
     $               alamb,iradia_diffra,xmed,ymed,anx,any,d,amu)
      implicit real*8(a-h,o-z)
	parameter (nmax=3000,nmax_eq=2*nmax)
	dimension fi(nmax),dfidn(nmax),
     $          xmed(nmax),ymed(nmax),anx(nmax),any(nmax)
      dimension AA(nmax,nmax),C(nmax,nmax,2),IPVT(nmax),b(nmax)

      DIMENSION Y(NMAX_EQ),DYDX(NMAX_EQ),YOUT(NMAX_EQ),
     >YT(NMAX_EQ),DYT(NMAX_EQ),DYM(NMAX_EQ)

      HH=H*0.5d0
      H6=H/6.d0
      XH=X+HH

      DO 11 I=1,N
        YT(I)=Y(I)+HH*DYDX(I)
  11  CONTINUE
      k_deriv=2
      t_sub=XH
      CALL DERIVS(XH,YT,DYT,AA,C,IPVT,
     $            ntot,n1,n2,n3,n4,n5,
     $            omega,amp_sway,amp_heave,amp_wavemaker,amp_wave,
     $            pigr,grav,fi,t_sub,
     $            alamb,iradia_diffra,xmed,ymed,anx,any,d,amu)
      DO 12 I=1,N
        YT(I)=Y(I)+HH*DYT(I)
  12  CONTINUE
      k_deriv=3
      t_sub=XH
      CALL DERIVS(XH,YT,DYM,AA,C,IPVT,
     $            ntot,n1,n2,n3,n4,n5,
     $            omega,amp_sway,amp_heave,amp_wavemaker,amp_wave,
     $            pigr,grav,fi,t_sub,
     $            alamb,iradia_diffra,xmed,ymed,anx,any,d,amu)
      DO 13 I=1,N
        YT(I)=Y(I)+H*DYM(I)
        DYM(I)=DYT(I)+DYM(I)
  13  CONTINUE
      k_deriv=4
      t_sub=X+H
      CALL DERIVS(X+H,YT,DYT,AA,C,IPVT,
     $            ntot,n1,n2,n3,n4,n5,
     $            omega,amp_sway,amp_heave,amp_wavemaker,amp_wave,
     $            pigr,grav,fi,t_sub,
     $            alamb,iradia_diffra,xmed,ymed,anx,any,d,amu)
      DO 14 I=1,N
        YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.d0*DYM(I))
  14  CONTINUE
      RETURN
      END
c*************************************************************************