C XKM       SOURCE    CB215821  26/08/24    21:18:55     12622          
      SUBROUTINE XKM(NCOURB,TRAC,XK,XM,W,F)
C
C ===============================================================
C CE SOUS-PROGRAMME EST APPELE PAR ENDOM1.
C IL CALCULE LES COEFF. K ET M DE LA FONCTION DECRIVANT LA COURBE
C DE TRACTION: Y=K*X**(1/M) OU X EST UNE DEFORM. ET Y UNE CONTR.
C
C ENTREES:
C -------
C NCOURB        = NBR. DE PTS. DE LA COURBE DE TRACTION
C TRAC(2*NCOURB)= COURBE DE TRACTION
C
C SORTIES:
C -------
C XK= COEFF. K
C XM= COEFF. M
C ================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION TRAC(*),B(2),C(2),W(*)
      DIMENSION A(2,2),F(NCOURB,2)
C
      DATA TOL,MAXITE/1.D-3,15/
C
C
C ==========================================
C CAS OU LA COURBE DE TRACTION N'A QUE 2 PTS
C ==========================================
      IF (NCOURB.EQ.2) THEN
        XK=(TRAC(3)-TRAC(1))/(TRAC(4)-TRAC(2))
        XM=1.D0
        GOTO 100
      ENDIF

C =======================================
C METHODE DES MOINDRES CARRES PONDERES
C ON LINEARISE: LOG(Y)=LOG(K)+1/M*LOG(X)
C C'EST DE LA FORME: YY=C1*F1(X)+C2*F2(X)
C ON RESOUD: A(2,2)C(2)=B(2)
C =======================================
      SIGY=TRAC(1)
      DO 10 I=1,NCOURB
        W(I)=(TRAC(2*I-1)-SIGY)*(TRAC(2*I-1)-SIGY)
   10 CONTINUE
C
      CALL ZERO(A,2,2)
      CALL ZDANUL(B,2)
C
C ...............................................................
C DS. LA SUITE, LE 1ER PT. DE LA COURBE EST OMIS CAR TRAC(1)=SIGY
C ET TRAC(2)=0 CE QUI POSE PB. POUR LE CALCUL DES LOGARITHMES
C ...............................................................
      DO 30 I=2,NCOURB
           A(1,1)=A(1,1)+W(I)
           A(1,2)=A(1,2)+W(I)*LOG(TRAC(2*I))
           A(2,2)=A(2,2)+W(I)*LOG(TRAC(2*I))*LOG(TRAC(2*I))
           B(1)=B(1)+W(I)*LOG(TRAC(2*I-1)-SIGY)
           B(2)=B(2)+W(I)*LOG(TRAC(2*I-1)-SIGY)*LOG(TRAC(2*I))
   30 CONTINUE
      DET=A(1,1)*A(2,2)-A(1,2)*A(1,2)
      C(1)=(B(1)*A(2,2)-B(2)*A(1,2))/DET
      C(2)=(B(2)*A(1,1)-B(1)*A(1,2))/DET
      XK=EXP(C(1))
      XM=1.D0/C(2)
C
C ----------------------------------------------------------------
C SI LA METHODE DE GAUSS-NEWTON NE CONVERGE PAS, ON RETIENT, COMME
C SOLUTIONS DU PB., XK ET XM CALCULES PAR LA METHODE PRECEDENTE,
C QU'ON STOCKE DS. XK0 ET XM0:
C ----------------------------------------------------------------
      XK0=XK
      XM0=XM
C
C =======================
C METHODE DE GAUSS-NEWTON
C =======================
      NITER=0
   40 DO 50 I=2,NCOURB
      F(I,1)=TRAC(2*I)**(1.D0/XM)
      F(I,2)=-(XK/XM**2.D0)*LOG(TRAC(2*I))*(TRAC(2*I)**(1.D0/XM))
   50 CONTINUE
C
      DO 102 I=1,2
       DO 101 J=1,2
        A(I,J)=0.D0
        DO 60 K=2,NCOURB
         A(I,J)=A(I,J)+F(K,I)*F(K,J)
   60 CONTINUE
  101  CONTINUE
  102 CONTINUE
C
      DO 103 I=1,2
        B(I)=0.D0
        DO 70 K=2,NCOURB
          B(I)=B(I)+F(K,I)*(TRAC(2*K-1)-SIGY-XK*TRAC(2*K)**(1.D0/XM))
   70   CONTINUE
  103 CONTINUE
C
      DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      DXK=(B(1)*A(2,2)-B(2)*A(1,2))/DET
      DXM=(B(2)*A(1,1)-B(1)*A(2,1))/DET
C
C -------------------------------------------------------------
C CALCUL DE L'ERREUR AV. (Z1) ET APRES (Z2) LA MAJ. DE XK ET XM
C -------------------------------------------------------------
      Z1=0.D0
      DO 80 K=2,NCOURB
        ZZ1=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM))
        ZZ1=ZZ1*ZZ1
        Z1=Z1+ZZ1
   80 CONTINUE
C
      XK=XK+DXK
      XM=XM+DXM
C
      Z2=0.D0
      DO 90 K=2,NCOURB
        ZZ2=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM))
        ZZ2=ZZ2*ZZ2
        Z2=Z2+ZZ2
   90 CONTINUE
C
C ----------------------------
C TESTS D'ARRET DES ITERATIONS
C ----------------------------
      NITER=NITER+1
C      ARRET=ABS(Z2-Z1)/Z2
C      IF (ARRET.LT.TOL) GOTO 100
      ARRET=ABS(Z2-Z1)
      IF (ARRET.LT.(TOL*Z2)) GOTO 100
      IF (NITER.LT.MAXITE) GOTO 40
C
C ---------------------------------------------------------------
C SI ON PASSE ICI, C'EST QUE G-NEWTON N'A PAS CONVERGE AU BOUT DE
C MAXITE ITER.: LA SOL. EST CELLE CALCULEE PAR LA METHODE DES
C MOINDRES CARRES PONDERES
C ---------------------------------------------------------------
      XK=XK0
      XM=XM0
  100 CONTINUE
      RETURN
      END