PROGRAM test_4
IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp)::coord(4,2),km(8,8),ym,pr,start_time,end_time
 INTEGER::i,j
 OPEN(10,FILE='fe95.dat')
 OPEN(11,FILE='fe95.res')
 WRITE(11,'(a/)')"      PLANE STRAIN 4-NODE QUADRILATERAL"
 WRITE(11,'(a)')"      Nodal Coordinates"
 DO i=1,4
   READ(10,*)coord(i,:)
   WRITE(11,'(2F12.4)')coord(i,:)
 END DO
 WRITE(11,'(/a)')"      Young's Modulus and Poisson's Ratio"
 READ(10,*)ym,pr
 WRITE(11,'(/2F12.4)')ym,pr
 CALL cpu_time(start_time)
! DO j=1,10000000
   CALL km_4node(coord,ym,pr,km)
! END DO
 WRITE(11,'(/a)')"      Stiffness Matrix"
 DO i=1,8
 DO j=1,i
   WRITE(11,'(2I3,F18.4)')i,j,km(i,j)
 END DO
 END DO
 CALL cpu_time(end_time)
 WRITE(11,'(/,A,F12.4)')"time taken=",end_time-start_time
STOP
CONTAINS

SUBROUTINE km_4node(coord,ym,pr,km)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::coord(4,2),ym,pr
 REAL(iwp),INTENT(OUT)::km(8,8)
 REAL(iwp)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp)::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec
 REAL(iwp)::f1,f2,f3,sm,ec1,ec2
 REAL(iwp)::d0=0.0_iwp,d1=1.0_iwp,d2=2.0_iwp

 x1=coord(1,1)
 y1=coord(1,2)
 x2=coord(2,1)
 y2=coord(2,2)
 x3=coord(3,1)
 y3=coord(3,2)
 x4=coord(4,1)
 y4=coord(4,2)

 sm=ym/d2/(d1+pr)
 ec1=ym*(d1-pr)/(d1+pr)/(d1-d2*pr)
 ec2=ec1*pr/(d1-pr)

!Term k11

 CALL family_a(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,g1,g2,   &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(1,1))

!Term k21

 CALL family_b(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)


 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(2,1))

!Term k22

 CALL family_a(y1,y4,y3,y2,x1,x4,x3,x2,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(2,2))

!Term k31

 CALL family_c(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(3,1))

!Term k32

 CALL family_d(y1,y2,y3,y4,x1,x2,x3,x4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL group_d(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(3,2))

!Term k33

 CALL family_a(x2,x3,x4,x1,y2,y3,y4,y1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(3,3))

!Term k41

 CALL family_d(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(4,1))

!Term k42

 CALL family_c(y2,y1,y4,y3,x2,x1,x4,x3,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(4,2))

!Term k43

 CALL family_b(x2,x3,x4,x1,y2,y3,y4,y1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(4,3))

!Term k44

 CALL family_a(y2,y1,y4,y3,x2,x1,x4,x3,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(4,4))

!Term k51

 CALL family_e(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(5,1))

!Term k52

 CALL family_f(y1,y2,y3,y4,x1,x2,x3,x4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL group_d(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(5,2))

!Term k53

 CALL family_c(x2,x3,x4,x1,y2,y3,y4,y1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(5,3))

!Term k54

 CALL family_d(y2,y3,y4,y1,x2,x3,x4,x1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL group_d(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(5,4))

!Term k55

 CALL family_a(x3,x4,x1,x2,y3,y4,y1,y2,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(5,5))

!Term k61

 CALL family_f(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(6,1))

!Term k62

 CALL family_e(y3,y2,y1,y4,x3,x2,x1,x4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(6,2))

!Term k63

 CALL family_d(x2,x3,x4,x1,y2,y3,y4,y1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(6,3))

!Term k64

 CALL family_c(y3,y2,y1,y4,x3,x2,x1,x4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(6,4))

!Term k65

 CALL family_b(x3,x4,x1,x2,y3,y4,y1,y2,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(6,5))

!Term k66

 CALL family_a(y3,y2,y1,y4,x3,x2,x1,x4,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(6,6))

!Term k71

 CALL family_c(x4,x1,x2,x3,y4,y1,y2,y3,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(7,1))

!Term k72

 CALL family_d(x4,x1,x2,x3,y4,y1,y2,y3,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(7,2))

!Term k73

 CALL family_e(x2,x3,x4,x1,y2,y3,y4,y1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(7,3))

!Term k74

 CALL family_f(y2,y3,y4,y1,x2,x3,x4,x1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL group_d(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(7,4))

!Term k75

 CALL family_c(x3,x4,x1,x2,y3,y4,y1,y2,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(7,5))

!Term k76

 CALL family_d(y3,y4,y1,y2,x3,x4,x1,x2,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL group_d(x3,x4,x1,x2,y3,y4,y1,y2,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(7,6))

!Term k77

 CALL family_a(x4,x1,x2,x3,y4,y1,y2,y3,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(7,7))

!Term k81

 CALL family_d(y4,y1,y2,y3,x4,x1,x2,x3,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL group_d(x4,x1,x2,x3,y4,y1,y2,y3,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(8,1))

!Term k82

 CALL family_c(y1,y4,y3,y2,x1,x4,x3,x2,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(8,2))

!Term k83

 CALL family_f(x2,x3,x4,x1,y2,y3,y4,y1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(8,3))

!Term k84

 CALL family_e(y4,y3,y2,y1,x4,x3,x2,x1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(8,4))

!Term k85

 CALL family_d(x3,x4,x1,x2,y3,y4,y1,y2,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(8,5))

!Term k86

 CALL family_c(y4,y3,y2,y1,x4,x3,x2,x1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(8,6))

!Term k87

 CALL family_b(x4,x1,x2,x3,y4,y1,y2,y3,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec2,f1,f2,f3,sm,km(8,7))

!Term k88

 CALL family_a(y4,y3,y2,y1,x4,x3,x2,x1,e1,e2,e3,e4,e5,e6,g1,g2,  &
   g3,g4,g5,g6,f1,f2,f3)

 CALL kterm (e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec1,f1,f2,f3,sm,km(8,8))

RETURN

END SUBROUTINE km_4node

SUBROUTINE kterm(e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec,f1,f2,f3,sm,k)
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN) ::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6,ec,f1,f2,f3,sm
 REAL(iwp),INTENT(OUT)::k
 REAL(iwp)::mm,pp,mp,pm,cons,l1m3,l1p3,l1m2,l1p2
 REAL(iwp)::dmm,dmp,dpp,dpm
 REAL(iwp)::d0=0.0_iwp,d1=1.0_iwp,d2=2.0_iwp,d3=3.0_iwp
 REAL(iwp)::d4=4.0_iwp,d6=6.0_iwp,d12=12.0_iwp,d48=48.0_iwp

 IF(f2.NE.d0.AND.f3.NE.d0)THEN

   mm=                         &
   -d2*ec*e3*e6*f2**5          &
   -d2*sm*g3*g6*f2**5          &
   -d2*ec*e2*e5*f3**5          &
   -d2*sm*g2*g5*f3**5          &
                                
   +d6*sm*g3*g6*f2**4*f1       &
   +d6*ec*e3*e6*f2**4*f1       &
   -d3*ec*e1*e6*f2**4*f3       &
   -d3*sm*g1*g6*f2**4*f3       &
   -d3*sm*g3*g4*f2**4*f3       &
   -d3*ec*e3*e4*f2**4*f3       &
   +d2*ec*e2*e6*f2**4*f3       &
   +d2*ec*e3*e5*f2**4*f3       &                         
   +d2*sm*g2*g6*f2**4*f3       &
   +d2*sm*g3*g5*f2**4*f3       &
   
   +d6*sm*g2*g5*f3**4*f1       &
   +d6*ec*e2*e5*f3**4*f1       &
   -d3*ec*e1*e5*f3**4*f2       &
   -d3*ec*e2*e4*f3**4*f2       &
   -d3*sm*g2*g4*f3**4*f2       &
   -d3*sm*g1*g5*f3**4*f2       &
   +d2*sm*g3*g5*f3**4*f2       &
   +d2*sm*g2*g6*f3**4*f2       &
   +d2*ec*e3*e5*f3**4*f2       &
   +d2*ec*e2*e6*f3**4*f2       &
                                
   +d2*sm*g3*g6*f1**3*f2**2    &
   +d2*ec*e3*e6*f1**3*f2**2    &
   +d2*sm*g2*g5*f1**3*f3**2    &
   +d2*ec*e2*e5*f1**3*f3**2    &
                                
   -d6*sm*g3*g6*f2**3*f1**2    &
   -d6*ec*e3*e6*f2**3*f1**2    &
   -d6*ec*e1*e4*f2**3*f3**2    &
   -d6*sm*g1*g4*f2**3*f3**2    &
   +d3*ec*e2*e4*f2**3*f3**2    &
   +d3*ec*e1*e5*f2**3*f3**2    &
   +d3*sm*g2*g4*f2**3*f3**2    &
   +d3*sm*g1*g5*f2**3*f3**2    &
   -d2*ec*e2*e5*f2**3*f3**2    &
   -d2*sm*g2*g5*f2**3*f3**2    &
                                
   -d6*ec*e2*e5*f3**3*f1**2    &
   -d6*sm*g2*g5*f3**3*f1**2    &
   -d6*ec*e1*e4*f3**3*f2**2    &
   -d6*sm*g1*g4*f3**3*f2**2    &
   +d3*sm*g1*g6*f3**3*f2**2    &
   +d3*sm*g3*g4*f3**3*f2**2    &
   +d3*ec*e1*e6*f3**3*f2**2    &
   +d3*ec*e3*e4*f3**3*f2**2    &
   -d2*ec*e3*e6*f3**3*f2**2    &
   -d2*sm*g3*g6*f3**3*f2**2    &
                                
   +sm*g2*g6*f1**3*f2*f3       &
   +sm*g3*g5*f1**3*f2*f3       &
   +ec*e2*e6*f1**3*f2*f3       &
   +ec*e3*e5*f1**3*f2*f3       &
                                
   +d6*ec*e3*e4*f2**3*f3*f1    &
   +d6*ec*e1*e6*f2**3*f3*f1    &
   +d6*sm*g3*g4*f2**3*f3*f1    &
   +d6*sm*g1*g6*f2**3*f3*f1    &
   -d3*ec*e3*e5*f2**3*f3*f1    &
   -d3*ec*e2*e6*f2**3*f3*f1    &
   -d3*sm*g2*g6*f2**3*f3*f1    &
   -d3*sm*g3*g5*f2**3*f3*f1    &
                                
   +d6*ec*e2*e4*f3**3*f1*f2    &
   +d6*sm*g2*g4*f3**3*f1*f2    &
   +d6*sm*g1*g5*f3**3*f2*f1    &
   +d6*ec*e1*e5*f3**3*f2*f1    &
   -d3*sm*g3*g5*f3**3*f1*f2    &
   -d3*ec*e2*e6*f3**3*f1*f2    &
   -d3*sm*g2*g6*f3**3*f1*f2    &
   -d3*ec*e3*e5*f3**3*f2*f1    &

   +d6*ec*e1*e4*f2**2*f3**2*f1  &
   +d6*sm*g1*g4*f2**2*f3**2*f1  &
   -d3*sm*g3*g4*f1**2*f2**2*f3  &
   -d3*ec*e3*e4*f1**2*f2**2*f3  &
   -d3*ec*e1*e6*f1**2*f2**2*f3  &
   -d3*sm*g1*g6*f1**2*f2**2*f3  &
   -d3*sm*g2*g4*f3**2*f1**2*f2  &
   -d3*ec*e1*e5*f3**2*f1**2*f2  &
   -d3*ec*e2*e4*f3**2*f1**2*f2  &
   -d3*sm*g1*g5*f3**2*f1**2*f2

!   mm=-d2*f2**5*ec*e3*e6-d2*f3**5*ec*e2*e5+sm*g2*g6*f1**3*f2*f3-d3*f3**4 &                              
!   *ec*e1*e5*f2+d6*f1*f3**4*sm*g2*g5+d6*f1*ec*e2*e4*f2*f3**3-d6*f1**2*ec*e2*e5*  &                      
!   f3**3+d2*f3**4*ec*e2*e6*f2-d3*f1*sm*g3*g5*f2*f3**3-d6*f1**2*sm*g2*g5*f3**3-d3* &                     
!   f1*ec*e3*e5*f2**3*f3+d2*f2**4*ec*e3*e5*f3-d3*sm*g2*g4*f1**2*f2*f3**2-d3*sm*g3 &                      
!   *g4*f1**2*f2**2*f3-d3*ec*e3*e4*f1**2*f2**2*f3+d2*f1**3*sm*g2*g5*f3**2+d2*     &                      
!   f2**4*ec*e2*e6*f3-d3*ec*e1*e5*f1**2*f2*f3**2-d3*f1*ec*e2*e6*f2**3*f3+d6*      &                      
!   f2**4*sm*g3*g6*f1+d6*ec*e3*e4*f1*f2**3*f3+d3*ec*e2*e4*f2**3*f3**2+d6*sm*g1*g4 &                      
!   *f1*f2**2*f3**2-d3*f1**2*ec*e2*e4*f2*f3**2+d6*ec*e1*e6*f1*f2**3*f3+d3*ec*e1   &                      
!   *e5*f2**3*f3**2+d2*f3**4*sm*g3*g5*f2-d3*f1*ec*e2*e6*f2*f3**3-d6*ec*e1*e4*     &
!   f2**2*f3**3+d3*sm*g1*g6*f2**2*f3**3-d2*ec*e2*e5*f2**3*f3**2+d6*sm*g2*g4*f1*f2 &
!   *f3**3-d3*f1*sm*g2*g6*f2**3*f3+d3*sm*g3*g4*f2**2*f3**3-d2*f3**5*sm*g2*g5-d2*   &
!   f2**5*sm*g3*g6-d6*ec*e1*e4*f2**3*f3**2+d6*sm*g3*g4*f1*f2**3*f3-d6*sm*g1*g4*   &
!   f2**3*f3**2-d3*f3**4*ec*e2*e4*f2+d6*sm*g1*g6*f1*f2**3*f3-d3*f2**4*ec*e1*e6*   &
!   f3-d3*f3**4*sm*g2*g4*f2-d2*ec*e3*e6*f2**2*f3**3-d6*sm*g3*g6*f1**2*f2**3+d2*sm* &
!   g3*g6*f1**3*f2**2+d2*ec*e3*e6*f1**3*f2**2-d6*sm*g1*g4*f2**2*f3**3-d3*ec*e1*   &
!   e6*f1**2*f2**2*f3+d2*f2**4*sm*g2*g6*f3-d3*f2**4*sm*g1*g6*f3+d2*f1**3*ec*e2*   &
!   e5*f3**2+d3*sm*g2*g4*f2**3*f3**2-d3*f2**4*sm*g3*g4*f3-d6*ec*e3*e6*f1**2*      &
!   f2**3+d6*ec*e1*e4*f1*f2**2*f3**2-d3*sm*g1*g6*f1**2*f2**2*f3+d6*f2**4*ec*e3    &
!   *e6*f1-d3*f2**4*ec*e3*e4*f3+d2*f2**4*sm*g3*g5*f3+d3*sm*g1*g5*f2**3*f3**2-d3*   &
!   f1*sm*g3*g5*f2**3*f3-d3*sm*g1*g5*f1**2*f2*f3**2+d3*ec*e1*e6*f2**2*f3**3-d2*sm &
!   *g2*g5*f2**3*f3**2+ec*e3*e5*f1**3*f2*f3-d3*f1*sm*g2*g6*f2*f3**3+d2*f3**4*sm  &
!   *g2*g6*f2+d3*ec*e3*e4*f2**2*f3**3-d3*f1*ec*e3*e5*f2*f3**3+d6*sm*g1*g5*f1*f2   &
!   *f3**3-d3*f3**4*sm*g1*g5*f2+d6*f1*ec*e1*e5*f2*f3**3-d2*sm*g3*g6*f2**2*f3**3+d2 &
!   *f3**4*ec*e3*e5*f2+f1**3*sm*g3*g5*f2*f3+ec*e2*e6*f1**3*f2*f3+d6*f1*f3**4*   &
!   ec*e2*e5

   mp=+d2*f2**5*ec*e3*e6-d2*f3**5*ec*e2*e5-sm*g2*g6*f1**3*f2*f3+d3*       &
   f3**4*ec*e1*e5*f2-d6*f1*f3**4*sm*g2*g5+d6*f1*ec*e2*e4*f2*f3**3-d6*f1**2*      &
   ec*e2*e5*f3**3+d2*f3**4*ec*e2*e6*f2+d3*f1*sm*g3*g5*f2*f3**3-d6*f1**2*sm*g2*   &
   g5*f3**3+d3*f1*ec*e3*e5*f2**3*f3-d2*f2**4*ec*e3*e5*f3+d3*sm*g2*g4*f1**2*f2    &
   *f3**2+d3*sm*g3*g4*f1**2*f2**2*f3+d3*ec*e3*e4*f1**2*f2**2*f3-d2*f1**3*sm*g2*  &
   g5*f3**2-d2*f2**4*ec*e2*e6*f3+d3*ec*e1*e5*f1**2*f2*f3**2+d3*f1*ec*e2*e6*      &
   f2**3*f3-d6*f2**4*sm*g3*g6*f1-d6*ec*e3*e4*f1*f2**3*f3-d3*ec*e2*e4*f2**3*      &
   f3**2-d6*sm*g1*g4*f1*f2**2*f3**2+d3*f1**2*ec*e2*e4*f2*f3**2-d6*ec*e1*e6*f1*   &
   f2**3*f3-d3*ec*e1*e5*f2**3*f3**2+d2*f3**4*sm*g3*g5*f2+d3*f1*ec*e2*e6*f2*      &
   f3**3-d6*ec*e1*e4*f2**2*f3**3-d3*sm*g1*g6*f2**2*f3**3+d2*ec*e2*e5*f2**3*      &
   f3**2+d6*sm*g2*g4*f1*f2*f3**3+d3*f1*sm*g2*g6*f2**3*f3-d3*sm*g3*g4*f2**2*      &
   f3**3-d2*f3**5*sm*g2*g5+d2*f2**5*sm*g3*g6+d6*ec*e1*e4*f2**3*f3**2-d6*sm*g3*g4  &
   *f1*f2**3*f3+d6*sm*g1*g4*f2**3*f3**2+d3*f3**4*ec*e2*e4*f2-d6*sm*g1*g6*f1*     &
   f2**3*f3+d3*f2**4*ec*e1*e6*f3+d3*f3**4*sm*g2*g4*f2-d2*ec*e3*e6*f2**2*         &
   f3**3+d6*sm*g3*g6*f1**2*f2**3-d2*sm*g3*g6*f1**3*f2**2-d2*ec*e3*e6*f1**3*      &
   f2**2-d6*sm*g1*g4*f2**2*f3**3+d3*ec*e1*e6*f1**2*f2**2*f3-d2*f2**4*sm*g2*g6*   &
   f3+d3*f2**4*sm*g1*g6*f3-d2*f1**3*ec*e2*e5*f3**2-d3*sm*g2*g4*f2**3*f3**2+d3*    &
   f2**4*sm*g3*g4*f3+d6*ec*e3*e6*f1**2*f2**3-d6*ec*e1*e4*f1*f2**2*f3**2+d3*      &
   sm*g1*g6*f1**2*f2**2*f3-d6*f2**4*ec*e3*e6*f1+d3*f2**4*ec*e3*e4*f3-d2*f2**4*   &
   sm*g3*g5*f3-d3*sm*g1*g5*f2**3*f3**2+d3*f1*sm*g3*g5*f2**3*f3+d3*sm*g1*g5*f1**2*&
   f2*f3**2-d3*ec*e1*e6*f2**2*f3**3+d2*sm*g2*g5*f2**3*f3**2-ec*e3*e5*f1**3*f2*  &
   f3+d3*f1*sm*g2*g6*f2*f3**3+d2*f3**4*sm*g2*g6*f2-d3*ec*e3*e4*f2**2*f3**3+d3*f1* &
   ec*e3*e5*f2*f3**3+d6*sm*g1*g5*f1*f2*f3**3+d3*f3**4*sm*g1*g5*f2+d6*f1*ec*e1*e5*&
   f2*f3**3-d2*sm*g3*g6*f2**2*f3**3+d2*f3**4*ec*e3*e5*f2-f1**3*sm*g3*g5*f2*f3-ec&
   *e2*e6*f1**3*f2*f3-d6*f1*f3**4*ec*e2*e5

   pp=+d2*f2**5*ec*e3*e6+d2*f3**5*ec*e2*e5+sm*g2*g6*f1**3*f2*f3-d3*       &
   f3**4*ec*e1*e5*f2+d6*f1*f3**4*sm*g2*g5-d6*f1*ec*e2*e4*f2*f3**3+d6*f1**2*ec*   &
   e2*e5*f3**3-d2*f3**4*ec*e2*e6*f2-d3*f1*sm*g3*g5*f2*f3**3+d6*f1**2*sm*g2*g5*   &
   f3**3-d3*f1*ec*e3*e5*f2**3*f3-d2*f2**4*ec*e3*e5*f3-d3*sm*g2*g4*f1**2*f2*      &
   f3**2-d3*sm*g3*g4*f1**2*f2**2*f3-d3*ec*e3*e4*f1**2*f2**2*f3+d2*f1**3*sm*g2*   &
   g5*f3**2-d2*f2**4*ec*e2*e6*f3-d3*ec*e1*e5*f1**2*f2*f3**2-d3*f1*ec*e2*e6*      &
   f2**3*f3+d6*f2**4*sm*g3*g6*f1-d6*ec*e3*e4*f1*f2**3*f3+d3*ec*e2*e4*f2**3*      &
   f3**2+d6*sm*g1*g4*f1*f2**2*f3**2-d3*f1**2*ec*e2*e4*f2*f3**2-d6*ec*e1*e6*f1    &
   *f2**3*f3+d3*ec*e1*e5*f2**3*f3**2-d2*f3**4*sm*g3*g5*f2-d3*f1*ec*e2*e6*f2*     &
   f3**3+d6*ec*e1*e4*f2**2*f3**3+d3*sm*g1*g6*f2**2*f3**3+d2*ec*e2*e5*f2**3*      &
   f3**2-d6*sm*g2*g4*f1*f2*f3**3-d3*f1*sm*g2*g6*f2**3*f3+d3*sm*g3*g4*f2**2*      &
   f3**3+d2*f3**5*sm*g2*g5+d2*f2**5*sm*g3*g6+d6*ec*e1*e4*f2**3*f3**2-d6*sm*g3*g4  &
   *f1*f2**3*f3+d6*sm*g1*g4*f2**3*f3**2-d3*f3**4*ec*e2*e4*f2-d6*sm*g1*g6*f1*     &
   f2**3*f3-d3*f2**4*ec*e1*e6*f3-d3*f3**4*sm*g2*g4*f2+d2*ec*e3*e6*f2**2*         &
   f3**3+d6*sm*g3*g6*f1**2*f2**3+d2*sm*g3*g6*f1**3*f2**2+d2*ec*e3*e6*f1**3*      &
   f2**2+d6*sm*g1*g4*f2**2*f3**3-d3*ec*e1*e6*f1**2*f2**2*f3-d2*f2**4*sm*g2*      &
   g6*f3-d3*f2**4*sm*g1*g6*f3+d2*f1**3*ec*e2*e5*f3**2+d3*sm*g2*g4*f2**3*         &
   f3**2-d3*f2**4*sm*g3*g4*f3+d6*ec*e3*e6*f1**2*f2**3+d6*ec*e1*e4*f1*f2**2*      &
   f3**2-d3*sm*g1*g6*f1**2*f2**2*f3+d6*f2**4*ec*e3*e6*f1-d3*f2**4*ec*e3*e4*      &
   f3-d2*f2**4*sm*g3*g5*f3+d3*sm*g1*g5*f2**3*f3**2-d3*f1*sm*g3*g5*f2**3*f3-d3*sm* &
   g1*g5*f1**2*f2*f3**2+d3*ec*e1*e6*f2**2*f3**3+d2*sm*g2*g5*f2**3*f3**2+ec*     &
   e3*e5*f1**3*f2*f3-d3*f1*sm*g2*g6*f2*f3**3-d2*f3**4*sm*g2*g6*f2+d3*ec*e3*e4*   &
   f2**2*f3**3-d3*f1*ec*e3*e5*f2*f3**3-d6*sm*g1*g5*f1*f2*f3**3-d3*f3**4*sm*g1*   &
   g5*f2-d6*f1*ec*e1*e5*f2*f3**3+d2*sm*g3*g6*f2**2*f3**3-d2*f3**4*ec*e3*e5*      &
   f2+f1**3*sm*g3*g5*f2*f3+ec*e2*e6*f1**3*f2*f3+d6*f1*f3**4*ec*e2*e5

   pm=-d2*f2**5*ec*e3*e6+d2*f3**5*ec*e2*e5-sm*g2*g6*f1**3*f2*f3+d3*      &
   f3**4*ec*e1*e5*f2-d6*f1*f3**4*sm*g2*g5-d6*f1*ec*e2*e4*f2*f3**3+d6*f1**2       &
   *ec*e2*e5*f3**3-d2*f3**4*ec*e2*e6*f2+d3*f1*sm*g3*g5*f2*f3**3+d6*f1**2*sm*g2   &
   *g5*f3**3+d3*f1*ec*e3*e5*f2**3*f3+d2*f2**4*ec*e3*e5*f3+d3*sm*g2*g4*f1**2*     &
   f2*f3**2+d3*sm*g3*g4*f1**2*f2**2*f3+d3*ec*e3*e4*f1**2*f2**2*f3-d2*f1**3*sm*   &
   g2*g5*f3**2+d2*f2**4*ec*e2*e6*f3+d3*ec*e1*e5*f1**2*f2*f3**2+d3*f1*ec*e2*e6    &
   *f2**3*f3-d6*f2**4*sm*g3*g6*f1+d6*ec*e3*e4*f1*f2**3*f3-d3*ec*e2*e4*f2**3*     &
   f3**2-d6*sm*g1*g4*f1*f2**2*f3**2+d3*f1**2*ec*e2*e4*f2*f3**2+d6*ec*e1*e6*f1    &
   *f2**3*f3-d3*ec*e1*e5*f2**3*f3**2-d2*f3**4*sm*g3*g5*f2+d3*f1*ec*e2*e6*f2*     &
   f3**3+d6*ec*e1*e4*f2**2*f3**3-d3*sm*g1*g6*f2**2*f3**3-d2*ec*e2*e5*f2**3*      &
   f3**2-d6*sm*g2*g4*f1*f2*f3**3+d3*f1*sm*g2*g6*f2**3*f3-d3*sm*g3*g4*f2**2*      &
   f3**3+d2*f3**5*sm*g2*g5-d2*f2**5*sm*g3*g6-d6*ec*e1*e4*f2**3*f3**2+d6*sm*g3*g4  &
   *f1*f2**3*f3-d6*sm*g1*g4*f2**3*f3**2+d3*f3**4*ec*e2*e4*f2+d6*sm*g1*g6*f1*     &
   f2**3*f3+d3*f2**4*ec*e1*e6*f3+d3*f3**4*sm*g2*g4*f2+d2*ec*e3*e6*f2**2*f3**3-   &
   6*sm*g3*g6*f1**2*f2**3-d2*sm*g3*g6*f1**3*f2**2-d2*ec*e3*e6*f1**3*f2**2+d6*sm* &
   g1*g4*f2**2*f3**3+d3*ec*e1*e6*f1**2*f2**2*f3+d2*f2**4*sm*g2*g6*f3+d3*f2**4*   &
   sm*g1*g6*f3-d2*f1**3*ec*e2*e5*f3**2-d3*sm*g2*g4*f2**3*f3**2+d3*f2**4*sm*g3*g4 &
   *f3-d6*ec*e3*e6*f1**2*f2**3-d6*ec*e1*e4*f1*f2**2*f3**2+d3*sm*g1*g6*f1**2*     &
   f2**2*f3-d6*f2**4*ec*e3*e6*f1+d3*f2**4*ec*e3*e4*f3+d2*f2**4*sm*g3*g5*f3-d3*sm  &
   *g1*g5*f2**3*f3**2+d3*f1*sm*g3*g5*f2**3*f3+d3*sm*g1*g5*f1**2*f2*f3**2-d3*ec*  &
   e1*e6*f2**2*f3**3-d2*sm*g2*g5*f2**3*f3**2-ec*e3*e5*f1**3*f2*f3+d3*f1*sm*g2*  &
   g6*f2*f3**3-d2*f3**4*sm*g2*g6*f2-d3*ec*e3*e4*f2**2*f3**3+d3*f1*ec*e3*e5*f2*   &
   f3**3-d6*sm*g1*g5*f1*f2*f3**3+d3*f3**4*sm*g1*g5*f2-d6*f1*ec*e1*e5*f2*f3**3+d2* &
   sm*g3*g6*f2**2*f3**3-d2*f3**4*ec*e3*e5*f2-f1**3*sm*g3*g5*f2*f3-ec*e2*e6*    &
   f1**3*f2*f3-d6*f1*f3**4*ec*e2*e5

   cons=+d4*ec*e3*e5*f1*f2**2*f3**2+d4*sm*g2*g6*f1*f2**2*f3**2+12*ec*e3* &
   e4*f2**3*f3**2-16*f1*sm*g2*g5*f2*f3**3+12*ec*e2*e4*f2**2*f3**3+12*sm*g3*   &
   g4*f2**3*f3**2+d4*f1*sm*g3*g5*f2**2*f3**2-16*ec*e3*e6*f1*f2**3*f3+12*sm*g1  &
   *g5*f2**2*f3**3+d4*ec*e2*e6*f1*f2**2*f3**2+12*ec*e1*e5*f2**2*f3**3-16*f1    &
   *ec*e2*e5*f2*f3**3+12*sm*g2*g4*f2**2*f3**3+12*sm*g1*g6*f2**3*f3**2-16*sm*  &
   g3*g6*f1*f2**3*f3+12*ec*e1*e6*f2**3*f3**2

   dmm=f1-f2-f3
   dmp=f1-f2+f3
   dpp=f1+f2+f3
   dpm=f1+f2-f3

   k=(mm*LOG(dmm)+mp*LOG(dmp)+pp*LOG(dpp)+pm*LOG(dpm)+cons)/d48/f2**3/f3**3

 ELSE IF(f2==d0.AND.f3.NE.d0)THEN
   
   l1m3=LOG(f1-f3)
   l1p3=LOG(f1+f3)
   
   k=-d1/d12*(                 &
   -d3*l1m3*f1*ec*e1*e6*f3     &
   -d3*l1m3*f1*ec*e3*e4*f3     &
   -d3*l1m3*f1*sm*g1*g6*f3     &
   -d3*l1m3*f1*sm*g3*g4*f3     &
   +d3*l1m3*f1**2*sm*g3*g6     &
   +d3*l1m3*sm*g1*g4*f3**2     &
   +d3*l1m3*ec*e1*e4*f3**2     &
   +l1m3*sm*g2*g5*f3**2        &
   +l1m3*ec*e2*e5*f3**2        &
   +d6*f1*sm*g3*g6*f3          &
   -d6*sm*g3*g4*f3**2          &
   +d6*f1*ec*e3*e6*f3          &
   -d6*sm*g1*g6*f3**2          &
   -d6*ec*e1*e6*f3**2          &
   -d6*ec*e3*e4*f3**2          &
   +d3*l1m3*f1**2*ec*e3*e6     &
   -d3*l1p3*ec*e1*e4*f3**2     &
   +d3*l1p3*f1*ec*e3*e4*f3     &
   +d3*l1p3*f1*ec*e1*e6*f3     &
   -d3*l1p3*sm*g1*g4*f3**2     &
   -d3*l1p3*f1**2*sm*g3*g6     &
   +d3*l1p3*f1*sm*g3*g4*f3     &
   +d3*l1p3*f1*sm*g1*g6*f3     &
   -d3*l1p3*f1**2*ec*e3*e6     &
   -l1p3*sm*g2*g5*f3**2        &
   -l1p3*ec*e2*e5*f3**2        &
    )/f3**3 
   
 ELSE IF(f2.NE.d0.AND.f3==d0)THEN
 
   l1m2=LOG(f1-f2)
   l1p2=LOG(f1+f2)

   k=-d1/d12*(                 &
   -l1p2*sm*g3*g6*f2**2        &
   -l1p2*ec*e3*e6*f2**2        &
   +l1m2*ec*e3*e6*f2**2        &
   +l1m2*sm*g3*g6*f2**2        &
   +d3*l1p2*f1*ec*e2*e4*f2     &
   -d3*l1p2*sm*g1*g4*f2**2     &
   -d6*ec*e2*e4*f2**2          &
   +d3*l1p2*f1*ec*e1*e5*f2     &
   -d6*sm*g2*g4*f2**2          &
   +d3*l1p2*f1*sm*g1*g5*f2     &
   -d3*l1p2*f1**2*sm*g2*g5     &
   +d3*l1m2*sm*g1*g4*f2**2     &
   -d3*l1m2*f1*ec*e1*e5*f2     &
   +d3*l1m2*ec*e1*e4*f2**2     &
   -d3*l1m2*f1*sm*g1*g5*f2     &
   +d3*l1m2*f1**2*ec*e2*e5     &
   -d3*l1m2*f1*sm*g2*g4*f2     &
   +d6*f1*sm*g2*g5*f2          &
   -d3*l1m2*f1*ec*e2*e4*f2     &
   -d6*sm*g1*g5*f2**2          &
   -d3*l1p2*f1**2*ec*e2*e5     &
   -d6*ec*e1*e5*f2**2          &
   +d3*l1m2*f1**2*sm*g2*g5     &
   +d3*l1p2*f1*sm*g2*g4*f2     &
   -d3*l1p2*ec*e1*e4*f2**2     &
   +d6*f1*ec*e2*e5*f2          &
     )/f2**3 
   
 ELSE

   k=d1/d6*(                   &
   +ec*e3*e6+sm*g3*g6          &
   +ec*e2*e5+sm*g2*g5          &
   +d3*ec*e1*e4                &
   +d3*sm*g1*g4                &
   )/f1  

 END IF

RETURN
END SUBROUTINE kterm

SUBROUTINE family_a(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,   &
  g1,g2,g3,g4,g5,g6,f1,f2,f3)

 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)

 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6
 REAL(iwp),INTENT(OUT)::f1,f2,f3

 e1= y4-y2
 e2=-y4+y3
 e3=-y3+y2
 e4= y4-y2
 e5=-y4+y3
 e6=-y3+y2

 g1=-x4+x2
 g2=-x3+x4
 g3= x3-x2
 g4=-x4+x2
 g5=-x3+x4
 g6= x3-x2

 f1=(x2-x4)*(y1-y3)-(x1-x3)*(y2-y4)
 f2=(x1-x2)*(y3-y4)-(x3-x4)*(y1-y2)
 f3=(x2-x3)*(y4-y1)-(x4-x1)*(y2-y3)  
 
RETURN
END SUBROUTINE family_a

SUBROUTINE family_b(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,   &
  g1,g2,g3,g4,g5,g6,f1,f2,f3)

 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 
 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6
 REAL(iwp),INTENT(OUT)::f1,f2,f3
 
 e1= y4-y2
 e2=-y4+y3
 e3=-y3+y2
 e4=-x4+x2
 e5=-x3+x4
 e6= x3-x2
 
 g1=-x4+x2
 g2=-x3+x4
 g3= x3-x2
 g4= y4-y2
 g5=-y4+y3
 g6=-y3+y2
  
 f1=(x2-x4)*(y1-y3)-(x1-x3)*(y2-y4)
 f2=(x1-x2)*(y3-y4)-(x3-x4)*(y1-y2)
 f3=(x2-x3)*(y4-y1)-(x4-x1)*(y2-y3)

RETURN
END SUBROUTINE family_b

SUBROUTINE family_c(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,   &
  g1,g2,g3,g4,g5,g6,f1,f2,f3)
 
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)

 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4 
 REAL(iwp),INTENT(OUT)::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6
 REAL(iwp),INTENT(OUT)::f1,f2,f3

 e1= y4-y2
 e2=-y4+y3
 e3=-y3+y2
 e4=-y3+y1
 e5=-y3+y4
 e6=-y1+y4
 
 g1=-x4+x2
 g2=-x3+x4
 g3= x3-x2
 g4= x3-x1
 g5= x3-x4
 g6=-x4+x1
 
 f1=(x2-x4)*(y1-y3)-(x1-x3)*(y2-y4)
 f2=(x1-x2)*(y3-y4)-(x3-x4)*(y1-y2)
 f3=(x2-x3)*(y4-y1)-(x4-x1)*(y2-y3)

RETURN
END SUBROUTINE family_c

SUBROUTINE family_d(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,   &
  g1,g2,g3,g4,g5,g6,f1,f2,f3)

 IMPLICIT NONE

 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4 
 REAL(iwp),INTENT(OUT)::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6
 REAL(iwp),INTENT(OUT)::f1,f2,f3

 e1= y4-y2
 e2=-y4+y3
 e3=-y3+y2
 e4= x3-x1
 e5= x3-x4
 e6=-x4+x1
 
 g1=-x4+x2
 g2=-x3+x4
 g3= x3-x2
 g4=-y3+y1
 g5=-y3+y4
 g6=-y1+y4
 
 f1=(x2-x4)*(y1-y3)-(x1-x3)*(y2-y4)
 f2=(x1-x2)*(y3-y4)-(x3-x4)*(y1-y2)
 f3=(x2-x3)*(y4-y1)-(x4-x1)*(y2-y3)

RETURN
END SUBROUTINE family_d

SUBROUTINE family_e(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,   &
  g1,g2,g3,g4,g5,g6,f1,f2,f3)

 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)

 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4 
 REAL(iwp),INTENT(OUT)::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6
 REAL(iwp),INTENT(OUT)::f1,f2,f3

 e1= y4-y2
 e2=-y4+y3
 e3=-y3+y2
 e4=-y4+y2
 e5=-y1+y2
 e6= y1-y4
 
 g1=-x4+x2
 g2=-x3+x4
 g3= x3-x2
 g4=-x2+x4
 g5=-x2+x1
 g6=-x1+x4
 
 f1=(x2-x4)*(y1-y3)-(x1-x3)*(y2-y4)
 f2=(x1-x2)*(y3-y4)-(x3-x4)*(y1-y2)
 f3=(x2-x3)*(y4-y1)-(x4-x1)*(y2-y3)

RETURN
END SUBROUTINE family_e

SUBROUTINE family_f(x1,x2,x3,x4,y1,y2,y3,y4,e1,e2,e3,e4,e5,e6,   &
  g1,g2,g3,g4,g5,g6,f1,f2,f3)

 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)

 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4 
 REAL(iwp),INTENT(OUT)::e1,e2,e3,e4,e5,e6,g1,g2,g3,g4,g5,g6
 REAL(iwp),INTENT(OUT)::f1,f2,f3

 e1= y4-y2
 e2=-y4+y3
 e3=-y3+y2
 e4=-x2+x4
 e5=-x2+x1
 e6=-x1+x4
 
 g1=-x4+x2
 g2=-x3+x4
 g3= x3-x2
 g4=-y4+y2
 g5=-y1+y2
 g6= y1-y4
 
 f1=(x2-x4)*(y1-y3)-(x1-x3)*(y2-y4)
 f2=(x1-x2)*(y3-y4)-(x3-x4)*(y1-y2)
 f3=(x2-x3)*(y4-y1)-(x4-x1)*(y2-y3)

RETURN
END SUBROUTINE family_f

SUBROUTINE group_d(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3)

 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)

 REAL(iwp),INTENT(IN)::x1,x2,x3,x4,y1,y2,y3,y4
 REAL(iwp),INTENT(OUT)::f1,f2,f3

 f1=(x2-x4)*(y1-y3)-(x1-x3)*(y2-y4)
 f2=(x1-x2)*(y3-y4)-(x3-x4)*(y1-y2)
 f3=(x2-x3)*(y4-y1)-(x4-x1)*(y2-y3)

RETURN
END SUBROUTINE group_d

END PROGRAM test_4
