program dynLLMax
	  
        implicit REAL*8 (A-H,O-Z)
      parameter Nmax=200,nHmax=50
 	logical erreur
      namelist/parametre/aa,xLfil,xMs,psi,Hk,sigma,Hsurface,Htorsion,
     &xnu,alpha,Aex
	namelist/donnee/courant,nformeI,frequence,instants,npoint,npsurf,
     &ngrad,ngradfonc,rapamort,nharm,ll,xprecision,Mperiode,nconveMax,
     &nRK,Hzhys1,npHys,nHzdyn,ndomaine,nbreHys,codemag
      namelist/champ/Hzi,Hzf,pasHz,Hzizone2,Hzfzone2,paszone2,
     &iretour,pascourant,courantf,nfrequence,Ffinal,pasfreq
	namelist/bobine/Hz0,nfHz,dephas

      dimension Vn(nHmax),harmcosa(nHmax),harmsina(nHmax),
     &harmcosb(nHmax),harmsinb(nHmax),Vnc(nhmax),Vns(nhmax)
      dimension dUfia(Nmax+2),dUza(Nmax+2),dUra(Nmax+2),Ri(Nmax+2)
      dimension dUfib(Nmax+2),dUzb(Nmax+2),dUrb(Nmax+2)
      dimension Ufia(Nmax+2),Uza(Nmax+2),Ura(Nmax+2)
      dimension Ufib(Nmax+2),Uzb(Nmax+2),Urb(Nmax+2)
      dimension Her(Nmax+2),Hefi(Nmax+2),Hez(Nmax+2)
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur
      common/para1/Hk,Hz,psi,xnu,xmu0,pi,pi2
	common/eMax/frequence,sigma,xMs,aa,courant,nformeI
      common/paraRK/deltaT,nRK
      common/para2/ani(Nmax+2),Htorsion,Hzhys1,xprecision,Dex,
     &FMuA,xLfil,pourdoma,codemag,instants,nharm,npHys,
     &Mperiode,nconveMax,npsurf,nHzdyn,ndomaine,nbreHys
	common/Seb/Hz0,omegaHz,dephas


      open(unit=7,file='DSeddaoui_donnees.in',status='old')
 
      open(unit=9,file='Vnf.dat',status='unknown')
      open(unit=10,file='VvsT.dat',status='unknown')
      open(unit=11,file='MzvsRvsTa.dat',status='unknown')
      open(unit=12,file='Mzmoyen.dat',status='unknown')
      open(unit=13,file='MzvsRvsTb.dat',status='unknown')
      open(unit=14,file='MfivsRvsTa.dat',status='unknown')
      open(unit=15,file='rapport.dat',status='unknown')
      open(unit=16,file='MfivsRvsTb.dat',status='unknown')
      open(unit=17,file='HfivsRvsTa.dat',status='unknown')
      open(unit=20,file='HfivsRvsTb.dat',status='unknown')
      open(unit=21,file='harm1cos.dat',status='unknown')
      open(unit=22,file='harm1sin.dat',status='unknown')
      open(unit=23,file='harm2cos.dat',status='unknown')
      open(unit=24,file='harm2sin.dat',status='unknown')
      open(unit=25,file='harm3cos.dat',status='unknown')
      open(unit=26,file='harm3sin.dat',status='unknown')
      open(unit=27,file='HzvsRvsTa.dat',status='unknown')
      open(unit=28,file='HzvsRvsTb.dat',status='unknown')
      open(unit=29,file='EfivsRvsTa.dat',status='unknown')
      open(unit=30,file='EfivsRvsTb.dat',status='unknown')
      open(unit=8,file='VfivsT.dat',status='unknown')
      open(unit=31,file='VfivsH.dat',status='unknown')
      open(unit=32,file='EzvsRvsTa.dat',status='unknown')
      open(unit=33,file='EzvsRvsTb.dat',status='unknown')
      open(unit=34,file='deriveflux.dat',status='unknown')
    
	read(7,parametre)
	read(7,donnee)
	read(7,champ)
	read(7,bobine)
   
      write(15,*)
      write(15,*) '          MMM  MMM  II  GGGGGG  NN  NN  LL'
      write(15,*) '          MM MM MM  II  GG      NNN NN  LL'
	write(15,*) '          MM MM MM  II  GG GGG  NNNNNN  LL'
	write(15,*) '          MM    MM  II  GG  GG  NN NNN  LL'
	write(15,*) '          MM    MM  II  GGGGGG  NN  NN  LLLLLL'
	write(15,*)
	write(15,*) 'Parametres du calcul'
	write(15,*) '==================================================='
	write(15,*)
	write(15,*) 'frequence      ', frequence/1.0d6,' MHz'
	write(15,*) 'courant        ',courant,'Arms'
	write(15,*) 'champ anisotropie ',Hk,' A/m'
	write(15,*) 'angle psi         ',psi,'degres'
	write(15,*) 'nbre de points    ',npoint
	write(34,*) ' Hz   f   dfia/dt   dfib/dt'
	if(pasHz.eq.0) then
	  write(15,*) 'champ Hz pour mesure',Hzi
	else
	  write(15,*) 'champ Hz pour mesure',Hzhys1
      endif
	write(15,*) 'nbre de point surface ',npsurf
	write(15,*) 'champ a la surface ',Hsurface
	write(15,*) 'ngrad  ', ngrad
	write(15,*) 'rayon du fil', aa
	write(15,*) 
	if(nHzdyn.eq.1) then
	  write(15,*) 'champ Hz dynamique inclus' 
	else
        write(15,*) 'champ Hz dynamique non inclus' 
	endif
	if(ndomaine.eq.1) then
	  write(15,*) ' 1 seul domaine' 
	else
        write(15,*) ' 2 domaines' 
	endif

      if (pasHz.ne.0.0d0) then
	  write(21,*) ' Hz  Re(V1a)   Re(V1b)'
	  write(22,*) ' Hz  Im(V1a)   Im(V1b)'
	  write(23,*) ' Hz  Re(V2a)   Re(V2b)'
	  write(24,*) ' Hz  Im(V2a)   Im(V2b)'
	  write(25,*) ' Hz  Re(V3a)   Re(V3b)'
	  write(26,*) ' Hz  Im(V3a)   Im(V3b)'
      endif

	pi=3.141593d0
	pi2=2.0d0*pi
      Rini=0.0d0
	Rfinal=1.0d0
	acare=aa*aa
	delR=1.0d0/(npoint+1.0d0)
	h=delR
	xmu0=pi*4.0d-7
	psi=psi*pi/1.8d2
	cospsi=dcos(psi)
	sinpsi=dsin(psi)
	omega=pi2*frequence
      deltaT=1.0d0/instants
	Dex=2.0d0*Aex/(xmu0*xMs*Hk*acare)
     	const=-xnu*xmu0*Hk/frequence
      FMuA=frequence*xmu0*aa
      write(10,*) '  t/T    I(A)    V(v)'
	write(31,*) ' Hz     Vfi(mVrms)'
	write(8,*) 't/T  VfiA(mV) VfiB(mV) Vfi(mV)'
      nbonboo=0
      codemag=codemag/Hk
	pourdoma=0.5d0
	if(ndomaine.eq.1) pourdoma=1.0d0
      if (nformeI.eq.1) courant=courant*1.414214d0

      omegaHz=pi2*nfHz
	dephas=dephas*pi2/180.0d0

      np2=npoint+2
	np1=npoint+1
	npint=np2-npsurf
	npint1=npint+1

	do i=npint1,np2
	  Ura(i)=0.0d0
	  Urb(i)=0.0d0
	  Ufia(i)=-dcos(psi)
 	  Ufib(i)=dcos(psi)
	  Uza(i)=-dsin(psi)
	  Uzb(i)=dsin(psi)
      enddo

 	do i=1,npoint+2
	 Ri(i)=(i-1)*delR
	enddo

 5	direcini=0.0d0
      if(Hzi.le.0.0) direcini=1.0d0
      do i=2,npint
	 angle=(pi/2.0d0-psi)*Ri(i)+direcini*pi
	 Uza(i)=dcos(angle)
	 Ura(i)=0.0d0
	 Ufia(i)=dsin(angle)
	 dUza(i)=0.0d0
	 dUra(i)=0.0d0
	 dUfia(i)=0.0d0
	 Uzb(i)=dcos(angle)
	 Urb(i)=0.0d0
	 Ufib(i)=dsin(angle)
	 dUzb(i)=0.0d0
	 dUrb(i)=0.0d0
	 dUfib(i)=0.0d0
	enddo
	if(Hzi.le.0.0d0) then
	  Uza(1)=-1.0d0
	  Uzb(1)=-1.0d0
	else
	  Uza(1)=1.0d0
	  Uzb(1)=1.0d0
	endif
      if(pasHz.ne.0.0d0) write(9,400) (n,n=1,nharm)

c====================================================================================
c calcul du gradient de l'anisotropie de surface qui s'étant sur ngrad intervalles entre points.

      
	n1g=npint1-ngrad
	n2g=npint1
      
	do i=1,n1g
	  ani(i)=0.0d0
	enddo
	do i=n2g,np2
	  ani(i)=Hsurface
	enddo

	if(ngrad.eq.0) goto 3

      if(ngradfonc.ne.1) goto 4
	x1g=Ri(n1g)
	x2g=Ri(n2g)
	y1g=0.0d0
	y2g=Hsurface
	axg=2.0d0/(x2g-x1g)
      bxg=(x2g+x1g)/(x1g-x2g)
	ayg=2.0d0/(y2g-y1g)
	byg=(y2g+y1g)/(y1g-y2g)
	
	do i=n1g,n2g-1
	  xg=Ri(i)
	  terme1g=-byg-0.5*(axg*xg+bxg)**3.0d0
	  terme2g=1.5d0*(axg*xg+bxg)
	  ani(i)=(terme1g+terme2g)/ayg

      enddo
 
  4   if(ngradfonc.ne.2) goto 3

      do i=1,npint1
        xg=Ri(i)
	  x0g=Ri(npint1)
	  ani(i)=Hsurface*exp(ngrad*(xg-x0g))
	enddo
c===========================================================================

 3    nbalay=0
	abpas1=dabs(pasHz)/2.0d0
	abpas2=dabs(paszone2)/2.0d0
      if(Hzi.gt.Hzf) nbalay=1
c nbalay=0 --> balayage de Hz de - vers +
c nbalay=1 --> balayage de Hz de + vers -
	if(nbalay.eq.0) then
	  if(pasHz.lt.0.0d0) pasHz=-pasHz 
	  if(paszone2.lt.0.0d0) paszone=-paszone2
	  if(Hzizone2.gt.Hzfzone2) then
	    Htemp=Hzizone2
	    Hzizone2=Hzfzone2
	    Hzfzone2=Htemp
	  endif
	endif
	if(nbalay.eq.1) then
	  if(pasHz.gt.0.0d0) pasHz=-pasHz 
	  if(paszone2.gt.0.0d0) paszone2=-paszone2
	  if(Hzizone2.lt.Hzfzone2) then
	    Htemp=Hzizone2
	    Hzizone2=Hzfzone2
	    Hzfzone2=Htemp
	  endif
	endif
      if(pasHz.eq.0.0d0) Hzhys1=Hzi

  1	Hz=Hzi

      if(nfrequence.ne.0) goto 10
  2   izone2=0
      if(nbalay.eq.0) then
        if((Hz.ge.Hzizone2) .and. (Hz.lt.Hzfzone2)) izone2=1
	endif
      if(nbalay.eq.1) then
        if((Hz.le.Hzizone2) .and. (Hz.gt.Hzfzone2)) izone2=1
	endif


      call dynamic(Ri,Ufia,Uza,Ura,dUfia,dUza,dUra,harmcosa,
     &harmsina,Ufib,Uzb,Urb,dUfib,dUzb,dUrb,harmcosb,harmsinb,
     &Vn,Vnc,Vns,tensionRMS)

      if (pasHz.ne.0.0d0) then
	  write(9,200) Hz,tensionRMS,Vn(1),Vn(2),Vn(3),Vn(4)
	  write(21,200) Hz,harmcosa(1),harmcosb(1)
	  write(22,200) Hz,harmsina(1),harmsinb(1)
	  write(23,200) Hz,harmcosa(2),harmcosb(2)
	  write(24,200) Hz,harmsina(2),harmsinb(2)
	  write(25,200) Hz,harmcosa(3),harmcosb(3)
	  write(26,200) Hz,harmsina(3),harmsinb(3)
	  if(izone2.eq.1) then
	    Hz=Hz+paszone2
          diffHz=dabs(Hz-Hzhys1)
   	    if(diffHz.lt.abpas2) Hz=Hzhys1 

	    goto 2
	  endif
        Hz=Hz+pasHz
	  diffHz=dabs(Hz-Hzhys1)
   	  if(diffHz.lt.abpas1) Hz=Hzhys1 

	  if((nbalay.eq.0) .and. (Hz.lt.Hzf)) goto 2
	  if((nbalay.eq.1) .and. (Hz.gt.Hzf)) goto 2
      endif


      if(nformeI.eq.1) then
        if(pascourant.ne.0.0d0) then
	    courf=courantf*1.414214d0
	    if(courant.lt.courf) then
	      courant=courant+pascourant*1.414214d0
	      goto 1
	    endif
	  endif
	endif
 
      if(nfrequence.eq.0) goto 12
 10   open(unit=18,file='VnvsF.dat',status='unknown')
      open(unit=19,file='VrVxvsF.dat',status='unknown')

      write(18,*) ' f   V1   V2   etc'
	write(19,*) ' f   Vrms   Vcosa    Vsina  Vcosb    Vsinb'
      Fini=frequence
 11   omega=pi2*frequence
     	const=-xnu*xmu0*Hk/frequence
      FMuA=frequence*xmu0*aa
      call dynamic(Ri,Ufia,Uza,Ura,dUfia,dUza,dUra,harmcosa,
     &harmsina,Ufib,Uzb,Urb,dUfib,dUzb,dUrb,harmcosb,harmsinb,
     &Vn,Vnc,Vns,tensionRMS)
	write(18,300) frequence,(Vn(i),i=1,nharm)
	write(19,300) frequence,tensionRMS,harmcosa(1),harmsina(1),
     &harmcosb(1),harmsinb(1)
	write(6,500) frequence
	write(6,*)
      if(frequence.lt.Ffinal) then
	  if(nfrequence.eq.1) frequence=frequence+pasfreq
	  if(nfrequence.eq.2) frequence=frequence*pasfreq
	  fff=frequence/Fini
	  if(fff.gt.2.0d0) then
	    Fini=frequence
	    if(instants.ge.5000) instants=int(instants/2)
	    deltaT=1.0d0/instants

	  endif
	  goto 11
	endif

 12    continue



      if(iretour.eq.0) goto 13
	iretour=0
	Hztemp=Hzi
	Hzi=Hzf
	Hzf=Hztemp
	pasHz=-pasHz
	Hztemp=Hzizone2
	Hzizone2=Hzfzone2
	Hzfzone2=Hztemp
	paszone2=-paszone2


	goto 5


 13   continue
      write(6,*) 'programme termine'
 100  format(2x,'  Hfi       Ufi      Uz      Ur     dUfi/dt')
 200  format(2x,f8.3,6(3x,f20.6))
 300  format(2x,f16.1,20(2x,f10.7))
 400  format(2x,'Hz',2x,'Vrms',2(2x,'Vn',I1),10(2x,I2))
 500  format(12x,'frequence=',E15.2)
	end

ccc  Fin du programme principal


c######################################################################################
c
      subroutine dynamic(Ri,Ufia,Uza,Ura,dUfia,dUza,dUra,harmcosa,
     &harmsina,Ufib,Uzb,Urb,dUfib,dUzb,dUrb,harmcosb,harmsinb,
     &Vn,Vnc,Vns,tensionRMS)
c
c
c    subroutine principale du calcul dynamique
c
c######################################################################################
       implicit REAL*8 (A-H,O-Z)

      parameter Nmax=200,nHmax=50
  	logical erreur
    
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur
      common/para1/Hk,Hz,psi,xnu,xmu0,pi,pi2
	common/eMax/frequence,sigma,xMs,aa,courant,nformeI
      common/paraRK/deltaT,nRK
      common/para2/ani(Nmax+2),Htorsion,Hzhys1,xprecision,Dex,
     &FMuA,xLfil,pourdoma,codemag,instants,nharm,npHys,
     &Mperiode,nconveMax,npsurf,nHzdyn,ndomaine,nbreHys
	common/Seb/Hz0,omegaHz,dephas

      dimension Hfia(Nmax+2),Hfib(Nmax+2),Hfimoin1(Nmax+2)
	dimension Ufia(Nmax+2),Uza(Nmax+2),Ura(Nmax+2)
	dimension Ufib(Nmax+2),Uzb(Nmax+2),Urb(Nmax+2)
      dimension Vn(nHmax),dHfia(Nmax+2),dHfib(Nmax+2),dBfi(Nmax+2)
      dimension alphaB(Nmax+1),bettaB(Nmax+1),gammaB(Nmax+1)
      dimension alphaU(Nmax+1),bettaU(Nmax+1),gammaU(Nmax+1)

      dimension alphai(Nmax+1),bettai(Nmax+1),gammai(Nmax+1)
	dimension ettaB(Nmax+1),ettai(Nmax+1),ettaU(Nmax+1)
      dimension deltafi(Nmax+2),deltar(Nmax+2),deltaz(Nmax+2)
	dimension Hefia(Nmax+2),Heza(Nmax+2),Hera(Nmax+2)
      dimension harmcosa(nHmax),harmsina(nHmax)
	dimension Hefib(Nmax+2),Hezb(Nmax+2),Herb(Nmax+2)
      dimension harmcosb(nHmax),harmsinb(nHmax)
      dimension Ufi1a(Nmax+2),Uz1a(Nmax+2),Ur1a(Nmax+2)
      dimension dUfia(Nmax+2),dUza(Nmax+2),dUra(Nmax+2),Ri(Nmax+2)
      dimension Ufi1b(Nmax+2),Uz1b(Nmax+2),Ur1b(Nmax+2)
      dimension dUfib(Nmax+2),dUzb(Nmax+2),dUrb(Nmax+2)
	dimension derfi1(Nmax+2),derfi2(Nmax+2),derz1(Nmax+2),
     &derz2(Nmax+2),derr1(Nmax+2),derr2(Nmax+2),somB0r(Nmax+2),
     &somBB(Nmax+2)
	dimension alphaUz(Nmax+1),bettaUz(Nmax+1),gammaUz(Nmax+1),
     &ettaUz(Nmax+1),Efi(Nmax+2),Vnc(nhmax),Vns(nhmax)
	dimension Efia(Nmax+2),Efib(Nmax+2),Eza(Nmax+2),Ezb(Nmax+2)
	dimension alphaEfi(Nmax+1),bettaEfi(Nmax+1),gammaEfi(Nmax+1),
     &ettaEfi(Nmax+1),Hzdyna(Nmax+2),Hzdynb(Nmax+2),rdUz(Nmax+2)
      dimension rUz(Nmax+2),alphaHz(Nmax+1),bettaHz(Nmax+1),
     &gammaHz(Nmax+1),ettaHz(Nmax+1),Hzdyn1a(Nmax+2),Hzdyn1b(Nmax+2)
	dimension alpharUz(Nmax+1),bettarUz(Nmax+1),gammarUz(Nmax+1),
     &ettarUz(Nmax+1),derHz2(Nmax+2)
      epsilon0=8.854188d-12
	aa2=aa*aa
	omega=pi2*frequence
      np2=npoint+2
	np1=npoint+1
	npint=np2-npsurf
	npint1=npint+1
      nperiode=0
      deltaT=1.0d0/instants
	instsauve=instants
      ns=1
c----------------------------------------------------------

c sauvegarde des valeurs des harmoniques de U et de ses composantes pour une comparaison après une periode
        do i=1,npoint+2
	    Ufi1a(i)=Ufia(i)
	    Uz1a(i)=Uza(i)
	    Ur1a(i)=Ura(i)
          dHfia(i)=0.0d0
	    Ufi1b(i)=Ufib(i)
	    Uz1b(i)=Uzb(i)
	    Ur1b(i)=Urb(i)
          dHfib(i)=0.0d0
	    Hzdyna(i)=0.0d0
	    Hzdynb(i)=0.0d0
	    Hzdyn1a(i)=0.0d0
	    Hzdyn1b(i)=0.0d0
	  enddo
	Vfia=0.0d0
	Vfib=0.0d0
      pourdom1=pourdoma
      cospsi=dcos(psi)
	sinpsi=dsin(psi)
	coef=epsilon0/(xmu0*Hk*Hk)

      nconverge=0
	nsauve=0
      npoinhys=int(instants/nbreHys)
      if(npoinhys.eq.0) npoinhys=100

c initialisation des coefficients trigonometriques de Fourier de U
	 do nh=1,nharm
	  harmcosa(nh)=0.0d0
	  harmsina(nh)=0.0d0
	  harmcosb(nh)=0.0d0
	  harmsinb(nh)=0.0d0
	enddo
	tensionRMS=0.0d0
	VfiRMS=0.0d0
      Uzmconva=0.0d0
	Uzmconvb=0.0d0
	pourdomb=1.0d0-pourdoma
	VinducaRMS=0.0d0
	VpeauaRMS=0.0d0
      VinducbRMS=0.0d0
	VpeaubRMS=0.0d0



c début de la période de temps.
c
  1   t=0.0d0
      dfluxmoya=0.0d0
      dfluxmoyb=0.0d0
      ninstant=0
	nstant=0
      Uzmoya=0.0d0
	Uzmoyb=0.0d0
c-----------------------------------------------------------


   2  continue
      ninstant=ninstant+1
	if(ninstant.eq.5000) then
	  write(6,600) t
	  ninstant=0
	endif


	if(pourdoma.ge.1.0d-3) then

c============================================================================================
c    calcul de la distribution du champ Hfi dans le fil à l'aide de l'équation de Maxwell-Ohm
c    pour le domaine a

        do np=1,npoint+2
	   Hfimoin1(np)=Hfia(np)
	  enddo
    
        call MaxwellOhm(t,dUfia,Hfimoin1,Hfia)

	  do k=1,npoint+2
	    dHfia(k)=(Hfia(k)-Hfimoin1(k))/deltaT
	  enddo
        if(nconverge.ne.0) then

c====================================================================
c calcul de la contribution en tension du domaine a
	    do np=1,npoint+2
	      dBfi(np)=dUfia(np)*xMs+dHfia(np)*Hk
	    enddo 
	    call interpole(Ri,dBfi,alphaB,bettaB,gammaB,ettaB,0)

c  somB0r(i) est l'integrale de dB/dt de Ri(1) à Ri(i) 
          call primitive(Ri,alphaB,bettaB,gammaB,ettaB,somB0r)

          do np=1,npoint+2
	      somBB(np)=2.0d0*Ri(np)*somB0r(np)

	    enddo
	    npp2=npoint+2

	    call interpole(Ri,somBB,alphai,bettai,gammai,ettai,0)

	    call integrale(0.0d0,1.0d0,Ri,alphai,bettai,gammai,ettai,Ez0a)

	    do np=1,npoint+2
            Eza(np)=somB0r(np)-Ez0a
	    enddo
          Vinduca=somB0r(npp2)*FMuA*xLfil
	    Vpeaua=-Ez0a*FMuA*xLfil
	    dfluxmoya=dfluxmoya+Vinduca*deltaT
c champ electrique axial et tension aux bornes du domaine a est:
	    if(nformeI.eq.1) xcourant=courant*dsin(pi2*t)
	    if(nformeI.eq.2) then
	      if(t.lt.5.0d-1) then
	        xcourant=courant
	      else
	        xcourant=-courant
	      endif
	    endif
	    if(nformeI.eq.3) then
	      if(t.lt.5.0d-2) then
	        xcourant=courant
	      else
	        xcourant=0.0d0
	      endif
	    endif


          Edc=xcourant/(pi*sigma*aa2)
  	    do i=1,npoint+2
	      Eza(i)=Eza(i)*FMuA+Edc
	    enddo
  	    tensiona=Eza(np2)*xLfil
   

 	    do nh=1,nharm
            factcos=dcos(nh*t*pi2)
	      factsin=dsin(nh*t*pi2)
	      harmcosa(nh)=harmcosa(nh)+2.0d0*tensiona*factcos*deltaT
	      harmsina(nh)=harmsina(nh)+2.0d0*tensiona*factsin*deltaT
	    enddo
        endif

calcul du champ dynamique Hz de la partie a

        if((nconverge.ne.0).or.(nHzdyn.eq.1)) then
    
          call MaxwellHz(t,dUza,Hzdyn1a,Hzdyna)

  	    call interpole(Ri,Hzdyna,alphaHz,bettaHz,gammaHz,ettaHz,0)
 	    call derive(Ri,alphaHz,bettaHz,gammaHz,ettaHz,Efi,derHz2)
  	    do i=2,npoint+2
            Efia(i)=-Efi(i)*Hk/(sigma*aa)
	    enddo
  	    Vfia=Efia(np2)*pi2*aa 
	    do i=1,np2
	      Hzdyn1a(i)=Hzdyna(i)
	    enddo
    
c-----------------------------------------------------------------------------
        endif
calcul de Uzma (la moyenne de la composante axiale de U) et de Uzmoya (moyenne temporelle de Uzma) pour le domaine a

        if((nconverge.ne.0).or.(codemag.ne.0)) then
  	    do k=1,npoint+2
            rUz(k)=Ri(k)*Uza(k)
	    enddo
          Uzma=0.0d0
  	    call interpole(Ri,rUz,alpharUz,bettarUz,gammarUz,ettarUz,0)
	    call integrale(0.0d0,1.0d0,Ri,alpharUz,bettarUz,gammarUz,
     &ettarUz,Uzma)
	    Uzma=Uzma*2.0d0
          Uzmoya=Uzmoya+Uzma*deltaT
        endif
calcul des derivées de U dans le domaine a
	  call interpole(Ri,Ufia,alphaU,bettaU,gammaU,ettaU,ns)
	  call derive(Ri,alphaU,bettaU,gammaU,ettaU,derfi1,derfi2)

 	  call interpole(Ri,Uza,alphaU,bettaU,gammaU,ettaU,ns)
	  call derive(Ri,alphaU,bettaU,gammaU,ettaU,derz1,derz2)

	  call interpole(Ri,Ura,alphaU,bettaU,gammaU,ettaU,ns)
	  call derive(Ri,alphaU,bettaU,gammaU,ettaU,derr1,derr2)



calcul des composantes du laplacien vecteur delta de U dans le domaine a.
        do i=2,npoint+2
	    rr=Ri(i)
	    rr2=rr*rr
          deltafi(i)=derfi2(i)+derfi1(i)/rr-Ufia(i)/rr2
          deltaz(i)=derz2(i)+derz1(i)/rr
	    deltar(i)=derr2(i)+derr1(i)/rr-Ura(i)/rr2
	  enddo
      

c  calcul du champ effectif He incluant Hz(appliqué), Hfi(sinusoidal circonférentiel), Ha(anisotropie), Hex(echange),
c     Htorsion et Hsurface dans le domaine a

        do i=2,npoint+2
          UscalNk=Ufia(i)*cospsi+Uza(i)*sinpsi
          UscalNtor=(Ufia(i)+Uza(i))*0.70710678
	      
	    Hefia(i)=Hfia(i)+UscalNk*cospsi+Dex*deltafi(i)
	    Heza(i)=Hz+UscalNk*sinpsi+Dex*deltaz(i)
	    Hera(i)=Dex*deltar(i)-xMs*Ura(i)/Hk

        enddo
	  if(nHzdyn.eq.1) then
	    do i=2,npoint+2
	      Heza(i)=Heza(i)+Hzdyna(i)
	    enddo
	  endif
        do i=2,npoint+2

	    Hefia(i)=Hefia(i)+5.0d-1*Htorsion*Ri(i)*(Ufia(i)+Uza(i))
	    Heza(i)=Heza(i)+5.0d-1*Htorsion*Ri(i)*(Ufia(i)+Uza(i))


	  enddo
       
	  do i=1,np2
          UscalNk=Ufia(i)*cospsi+Uza(i)*sinpsi
          Hefia(i)=Hefia(i)+ani(i)*UscalNk*cospsi
          Heza(i)=Heza(i)+ani(i)*UscalNk*sinpsi
        enddo
      endif







	if(pourdomb.ge.1.0d-3) then

c    calcul de la distribution du champ Hfi dans le fil à l'aide de l'équation de Maxwell-Ohm
c    pour le domaine b

        do np=1,npoint+2
	   Hfimoin1(np)=Hfib(np)
	  enddo
    
        call MaxwellOhm(t,dUfib,Hfimoin1,Hfib)

	  do k=1,npoint+2
	    dHfib(k)=(Hfib(k)-Hfimoin1(k))/deltaT
	  enddo




        if(nconverge.ne.0) then


c============================================================================
c calcul de la contribution en tension du domaine b
	    do np=1,npoint+2
	      dBfi(np)=dUfib(np)*xMs+dHfib(np)*Hk
	    enddo 
	    call interpole(Ri,dBfi,alphaB,bettaB,gammaB,ettaB,0)

c  somB0r(i) est l'integrale de dB/dt de Ri(1) à Ri(i) 
          call primitive(Ri,alphaB,bettaB,gammaB,ettaB,somB0r)

          do np=1,npoint+2
	      somBB(np)=2.0d0*Ri(np)*somB0r(np)
	    enddo
	    npp2=npoint+2

	    call interpole(Ri,somBB,alphai,bettai,gammai,ettai,0)
	    call integrale(0.0d0,1.0d0,Ri,alphai,bettai,gammai,ettai,Ez0b)

	    do np=1,npoint+2
            Ezb(np)=somB0r(np)-Ez0b
	    enddo
          Vinducb=somB0r(npp2)*FMuA*xLfil
	    Vpeaub=-Ez0b*FMuA*xLfil
	    dfluxmoyb=dfluxmoyb+Vinducb

c Champ electrique axial et tension aux bornes du domaine b:
	    if(nformeI.eq.1) xcourant=courant*dsin(pi2*t)
	    if(nformeI.eq.2) then
	      if(t.lt.5.0d-1) then
	        xcourant=courant
	      else
	        xcourant=-courant
	      endif
	    endif
	    if(nformeI.eq.3) then
	      if(t.lt.5.0d-2) then
	        xcourant=courant
	      else
	        xcourant=0.0d0
	      endif
	    endif

	    Edc=xcourant/(pi*sigma*aa2)
	    do i=1,npoint+2
	      Ezb(i)=Ezb(i)*FMuA+Edc
	    enddo
          tensionb=Ezb(np2)*xLFil

        
 	    do nh=1,nharm
            factcos=dcos(nh*t*pi2)
	      factsin=dsin(nh*t*pi2)
	      harmcosb(nh)=harmcosb(nh)+2.0d0*tensionb*factcos*deltaT
	      harmsinb(nh)=harmsinb(nh)+2.0d0*tensionb*factsin*deltaT
	    enddo
	  endif


calcul du champ dynamique Hz de la partie b

        if((nconverge.ne.0).or.(nHzdyn.eq.1)) then
 
          call MaxwellHz(t,dUzb,Hzdyn1b,Hzdynb)

  	    call interpole(Ri,Hzdynb,alphaHz,bettaHz,gammaHz,ettaHz,0)
 	    call derive(Ri,alphaHz,bettaHz,gammaHz,ettaHz,Efi,derHz2)
  	    do i=2,npoint+2
            Efib(i)=-Efi(i)*Hk/(sigma*aa)
	    enddo
   	    Vfib=Efib(np2)*pi2*aa 

  	    do i=1,np2
	      Hzdyn1b(i)=Hzdynb(i)
	    enddo
        endif
calcul de Uzmb (la moyenne de la composante axiale de U) et de Uzmoyb (moyenne temporelle de Uzmb) pour le domaine b

        if((nconverge.ne.0).or.(codemag.ne.0)) then

 	    do k=1,npoint+2
            rUz(k)=Ri(k)*Uzb(k)
	    enddo

          Uzmb=0.0d0
  	    call interpole(Ri,rUz,alpharUz,bettarUz,gammarUz,ettarUz,0)
	    call integrale(0.0d0,1.0d0,Ri,alpharUz,bettarUz,gammarUz,
     &ettarUz,Uzmb)
	    Uzmb=Uzmb*2.0d0
          Uzmoyb=Uzmoyb+Uzmb*deltaT
        endif




calcul des derivées de U dans le domaine b
	  call interpole(Ri,Ufib,alphaU,bettaU,gammaU,ettaU,ns)
	  call derive(Ri,alphaU,bettaU,gammaU,ettaU,derfi1,derfi2)

 	  call interpole(Ri,Uzb,alphaU,bettaU,gammaU,ettaU,ns)
	  call derive(Ri,alphaU,bettaU,gammaU,ettaU,derz1,derz2)

	  call interpole(Ri,Urb,alphaU,bettaU,gammaU,ettaU,ns)
	  call derive(Ri,alphaU,bettaU,gammaU,ettaU,derr1,derr2)


calcul des composantes du laplacien vecteur delta de U dans le domaine b.
        do i=2,npoint+2
	    rr=Ri(i)
	    rr2=rr*rr
          deltafi(i)=derfi2(i)+derfi1(i)/rr-Ufib(i)/rr2
          deltaz(i)=derz2(i)+derz1(i)/rr
	    deltar(i)=derr2(i)+derr1(i)/rr-Urb(i)/rr2
	  enddo
c        deltafi(np2)=(derfi1(np2)+Ufib(np2))*aa
c        deltaz(np2)=(derz1(np2)+Uzb(np2))*aa
c        deltar(np2)=(derr1(np2)+Urb(np2))*aa
      

c  calcul du champ effectif He incluant Hz(appliqué), Hfi(sinusoidal circonférentiel), Ha(anisotropie), Hex(echange),
c     Htorsion et Hsurface dans le domaine b

        do i=2,npoint+2
          UscalNk=Ufib(i)*cospsi+Uzb(i)*sinpsi
          UscalNtor=(Ufib(i)+Uzb(i))*0.70710678
	      
	    Hefib(i)=Hfib(i)+UscalNk*cospsi+Dex*deltafi(i)
	    Hezb(i)=Hz+UscalNk*sinpsi+Dex*deltaz(i)
	    Herb(i)=Dex*deltar(i)-xMs*Urb(i)/Hk

	    Hefib(i)=Hefib(i)+5.0d-1*Ri(i)*Htorsion*(Ufib(i)+Uzb(i))
	    Hezb(i)=Hezb(i)+5.0d-1*Htorsion*Ri(i)*(Ufib(i)+Uzb(i))
	  enddo
	  if(nHzdyn.eq.1) then
	    do i=2,npoint+2
	      Hezb(i)=Hezb(i)+Hzdynb(i)
	    enddo
	  endif


	  do i=1,np2
          UscalNk=Ufib(i)*cospsi+Uzb(i)*sinpsi
          Hefib(i)=Hefib(i)+ani(i)*UscalNk*cospsi
          Hezb(i)=Hezb(i)+ani(i)*UscalNk*sinpsi
        enddo

      endif

      tensiontot=pourdoma*tensiona+pourdomb*tensionb
      Vfi=pourdoma*Vfia+pourdomb*Vfib
	if(nconverge.ne.0) then
	  tensionRMS=tensionRMS+(tensiontot**2)*deltaT
        VfiRMS=VfiRMS+(Vfi**2)*deltaT
        VinducaRMS=VinducaRMS+(Vinduca**2)*deltaT
	  VpeauaRMS=VpeauaRMS+(Vpeaua**2)*deltaT
        VinducbRMS=VinducbRMS+(Vinducb**2)*deltaT
	  VpeaubRMS=VpeaubRMS+(Vpeaub**2)*deltaT

	endif


c contribusion du champ de désaimantation axial



      Uzm=pourdoma*Uzma+pourdomb*Uzmb
      do i=1,npoint+2
        Heza(i)=Heza(i)-codemag*Uzm
        Hezb(i)=Hezb(i)-codemag*Uzm
      enddo

conservation du flux traversant la paroie.(non pris en compte ici)
c      if((pourdoma.gt.1.0d-3).and.(pourdoma.lt.9.99d-1)) then
c	  pourdoma2=pourdoma*pourdoma
c	  pourdomb2=pourdomb*pourdomb
c        do i=1,npoint+2
c           Upera=Uza(i)*cospsi-Ufia(i)*sinpsi
c	     Uperb=Uzb(i)*cospsi-Ufib(i)*sinpsi
c           Heza(i)=Heza(i)+pourdomb2*xMs*(Uperb-Upera)*cospsi/Hk
c	     Hezb(i)=Hezb(i)+pourdoma2*xMs*(Upera-Uperb)*cospsi/Hk
c           Hefia(i)=Hefia(i)-pourdomb2*xMs*(Uperb-Upera)*sinpsi/Hk
c	     Hefib(i)=Hefib(i)-pourdoma2*xMs*(Upera-Uperb)*sinpsi/Hk
c        enddo
c	endif



	if(nconverge.ne.0) then
	  if(nstant.eq.npoinhys) then
    	    nstant=0
          if(Hz.eq.HzHys1) then

	      write(10,300) t,xcourant,tensiontot
	      write(8,900) t,Vfia,Vfib,Vfi
            write(11,100) (Uza(np),np=1,npoint+2)
            write(13,100) (Uzb(np),np=1,npoint+2)
            write(14,100) (Ufia(np),np=1,npoint+2)
            write(16,100) (Ufib(np),np=1,npoint+2)
            write(17,100) (Hfia(np),np=1,npoint+2)
            write(20,100) (Hfib(np),np=1,npoint+2)
            write(27,100) (Hzdyna(np),np=1,npoint+2)
            write(28,100) (Hzdynb(np),np=1,npoint+2)
            write(29,100) (Efia(np),np=1,npoint+2)
            write(30,100) (Efib(np),np=1,npoint+2)
            write(32,100) (Eza(np),np=1,npoint+2)
            write(33,100) (Ezb(np),np=1,npoint+2)

		endif

	  endif

      endif
  
c    calcul de dU/dt à l'instant t et la position de U à l'instant t+deltaT 
	inst=instants
	if(pourdoma.ge.1.0d-3) then
	  call RK4(t,inst,Ufia,Uza,Ura,Hefia,Heza,Hera,dUfia,dUza,dUra)
	endif
	if(pourdoma.le.0.999d0) then
	  call RK4(t,inst,Ufib,Uzb,Urb,Hefib,Hezb,Herb,dUfib,dUzb,dUrb)
      endif

	if(inst.ne.instants) then
	  instants=inst
        write(6,700) t,Hz,inst
        write(15,700) t,Hz,inst
	endif

      Uzab2=pourdoma*Uza(2)+pourdomb*Uzb(2)
      if(Uzab2.gt.0.1d0) then
	  Uza(1)=1.0d0
	  Uzb(1)=1.0d0
	endif
	if(Uzab2.lt.-0.1d0) then
	  Uza(1)=-1.0d0
	  Uzb(1)=-1.0d0
	endif


	t=t+deltaT

c   a ce stade, le temps vient d'étre incrémenté de deltaT
c   la nouvelle position de U(r) est obtenue pour t+deltaT



c  tant que la periode n'est pas complétée, on repart vers 2 pour determiner He et ainsi de suite
	nstant=nstant+1

	if(t.lt.1.0) goto 2

c à ce stade, la periode est bouclée.
      nperiode=nperiode+1


c reajustement du coeficien pourdom (pourcentage de la longeurs du domaine a par rapport à la longeur totale) (non pris en compte ici)

c      if(ndomaine.ne.1)
c        Hda=codemag*Uzmoya
c        Hdb=codemag*Uzmoyb
c	  difHd=dabs(Hda-Hdb)
c	  Hdaeps=1.001d0*Hda+1.0d-5
c	  Hdbeps=1.001d0*Hdb-1.0d-5
c 	  if(difHd.gt.1.0d-3) then
c	    difdom=(Hz-Hdb)/(Hda-Hdb)-pourdoma
c	    pourdoma=pourdoma+1.0d0*difdom
c	  endif
c        if((Hz.lt.Hdbeps).or.(pourdoma.le.1.0d-2)) pourdoma=0.0d0
c	  if((Hz.gt.Hdaeps).or.(pourdoma.ge.9.9d-1)) pourdoma=1.0d0
c        if(pourdoma.le.1.0d-3) pourdoma=0.0d0
c	  if(pourdoma.ge.9.99d-1) pourdoma=1.0d0

c	  pourdomb=1.0d0-pourdoma
c      endif



c  comparaison des composantes de U avant et après la periode
	somU=100.0d0*dabs(pourdoma-pourdom1)
      do i=1,npoint+2
	    somU=somU+dabs(Ufia(i)-Ufi1a(i))
	    somU=somU+dabs(Uza(i)-Uz1a(i))
	    somU=somU+dabs(Ura(i)-Ur1a(i))
	    somU=somU+dabs(Ufib(i)-Ufi1b(i))
	    somU=somU+dabs(Uzb(i)-Uz1b(i))
	    somU=somU+dabs(Urb(i)-Ur1b(i))
	enddo
      if(nconverge.ne.0) then
        Uzmconva=Uzmconva+Uzmoya
	  Uzmconvb=Uzmconvb+Uzmoyb
	  write(34,201) Hz,frequence,dfluxmoya,dfluxmoyb
	endif

	write(6,*)
      write(6,500) Hz,courant,somU
c	write(6,*) 'alpha = ',pourdoma
      write(6,*)

c si la convergence n'est pas atteinte (rapport > xpercision) refaire une autre periode
	if(nconverge.eq.nconveMax) goto 10
	if(nperiode.ge.Mperiode) then
	  nconverge=nconverge+1
	  write(15,400) Hz,courant,somU,pourdoma
	  if(somU.le.xprecision) then
	    nsauve=nconverge
	    nconverge=nconveMax
	  endif
	  goto1
	endif
      if(somU.gt.xprecision) then
        do i=1,npoint+2
	      Ufi1a(i)=Ufia(i)
	      Uz1a(i)=Uza(i)
	      Ur1a(i)=Ura(i)
	      Ufi1b(i)=Ufib(i)
	      Uz1b(i)=Uzb(i)
	      Ur1b(i)=Urb(i)
	  enddo
        pourdom1=pourdoma
	  goto 1
      endif
	 
c à ce stade, la convergence est atteinte. écriture des résultats.
      
      nconverge=nconveMax
  	goto 1
 10   continue
	if(nperiode.gt.Mperiode) then
	  if(nsauve.eq.0) nsauve=nconveMax
        do nh=1,nharm
          harmcosa(nh)=harmcosa(nh)/nsauve
          harmsina(nh)=harmsina(nh)/nsauve
          harmcosb(nh)=harmcosb(nh)/nsauve
          harmsinb(nh)=harmsinb(nh)/nsauve
	  enddo
	  tensionRMS=tensionRMS/nsauve
	  VfiRMS=VfiRMS/nsauve
	  VinducaRMS=VinducaRMS/nsauve
	  VpeauaRMS=VpeauaRMS/nsauve
        VinducbRMS=VinducbRMS/nsauve
	  VpeaubRMS=VpeaubRMS/nsauve

	  Uzmconva=Uzmconva/nsauve
	  Uzmconvb=Uzmconvb/nsauve
        write(15,*)
	endif
      
      do nh=1,nharm
	  Vnc(nh)=pourdoma*harmcosa(nh)+pourdomb*harmcosb(nh)
	  Vns(nh)=pourdoma*harmsina(nh)+pourdomb*harmsinb(nh)
	  Vn(nh)=dsqrt(Vnc(nh)**2+Vns(nh)**2)
	enddo

     	instants=instsauve

      tensionRMS=sqrt(tensionRMS)
	VfiRMS=sqrt(VfiRMS)


      write(12,300) Hz,Uzmconva,Uzmconvb,pourdoma
	write(31,800) Hz,VfiRMS

 100  format(102(3x,f10.5))
 200  format(2x,f12.6,4(3x,f10.5))
 201  format(2x,f12.6,2x,f15.1,2(2x,f10.5))
 300  format(2x,f12.6,4(3x,f10.5))
 400  format(' convergence non atteinte pour Hz=',f8.4,'courant=',f8.6,
     &'decalage=',f8.5,'portion du domaaine a =',f8.6)
 500  format('  Hz=',f10.3,'  courant(Ampl)=',f7.4,'  decalage=',f11.7)
 600  format('  temps= ',f7.5)
 700  format(' pas du temps change a t=',f8.6,' et Hz =',f8.4,'
     &pour la valeur ',i12) 
 800  format(2x,f12.6,3x,f12.6)
 900  format(2x,f12.6,3(3x,f12.8))
      return
 	end
	








      subroutine eqLLG(Ufi,Uz,Ur,Hefi,Hez,Her,dUfi,dUz,dUr)
 
c
c     subroutine de calcul de dU/dt selon l'equation de : Landau-Lifshitz si ll=1
c                                                         Gilbert         si ll#1
c
c  Ufi,Uz,Ur sont les composantes de U(r)=M(r)/Ms
c  Hefi,Hez,Her sont les composantes du champ effectif
c dUfi,dUz,dUr sont les composantes de dU(r)/dt
       implicit REAL*8 (A-H,O-Z)
      parameter Nmax=200
	logical erreur

      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur
      dimension Ufi(Nmax+2),Uz(Nmax+2),Ur(Nmax+2)
      dimension Hefi(Nmax+2),Hez(Nmax+2),Her(Nmax+2)
      dimension dUfi(Nmax+2),dUz(Nmax+2),dUr(Nmax+2)
      dimension UHfi(Nmax+2),UHz(Nmax+2),UHr(Nmax+2)
      dimension UUHfi(Nmax+2),UUHz(Nmax+2),UUHr(Nmax+2)


	call prodvect(Ufi,Uz,Ur,Hefi,Hez,Her,UHfi,UHz,UHr)
      
	np=npoint+2
	call prodvect(Ufi,Uz,Ur,UHfi,UHz,UHr,UUHfi,UUHz,UUHr)
	if(ll.eq.1) then
	  do n=1,npoint+1
 	    dUfi(n)=const*(UHfi(n)+alpha*UUHfi(n))
	    dUz(n)=const*(UHz(n)+alpha*UUHz(n))
	    dUr(n)=const*(UHr(n)+alpha*UUHr(n))
        enddo
	  constsurf=const/rapamort
 	    dUfi(np)=constsurf*UHfi(np)+const*alpha*UUHfi(np)
	    dUz(np)=constsurf*UHz(np)+const*alpha*UUHz(np)
	    dUr(np)=constsurf*UHr(np)+const*alpha*UUHr(np)

      else
        alphasurf=alpha*rapamort
	  alpha2=alpha*alpha
	  alpha2sur=alphasurf*alphasurf
	  constG=const/(1.0d0+alpha2)
	  constGsur=const/(1.0d0+alpha2sur)
	  do n=1,npoint+1
 	    dUfi(n)=constG*(UHfi(n)+alpha*UUHfi(n))
	    dUz(n)=constG*(UHz(n)+alpha*UUHz(n))
	    dUr(n)=constG*(UHr(n)+alpha*UUHr(n))
        enddo
	np=npoint+2
 	    dUfi(np)=constGsur*(UHfi(np)+alphasurf*UUHfi(np))
	    dUz(np)=constGsur*(UHz(np)+alphasurf*UUHz(np))
	    dUr(np)=constGsur*(UHr(np)+alphasurf*UUHr(np))

	endif
      return
	end




c===============================================================
c
c
      subroutine prodvect(X1,Y1,Z1,X2,Y2,Z2,XX,YY,ZZ)
c
c     subroutine produit vectoriel (XX,YY,ZZ)=(X1,Y1,Z1)x (X2,Y2,Z2)
c     où l'argument correspond au numero du point sur l'axe radial du fil
      implicit REAL*8 (A-H,O-Z)

      parameter Nmax=200
	logical erreur

	dimension X1(Nmax+2),Y1(Nmax+2),Z1(Nmax+2)
	dimension X2(Nmax+2),Y2(Nmax+2),Z2(Nmax+2)
	dimension XX(Nmax+2),YY(Nmax+2),ZZ(Nmax+2)
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur

	do n=1,npoint+2
	 XX(n)=Y1(n)*Z2(n)-Y2(n)*Z1(n)
	 YY(n)=X2(n)*Z1(n)-X1(n)*Z2(n)
	 ZZ(n)=X1(n)*Y2(n)-X2(n)*Y1(n)
	enddo
	return
	end







      subroutine MaxwellOhm(t,dUfi,Hi1,Hi)


c#################################################################################################################
c#                                                                                                               #
c#                 subroutine de calcul du champ magnétique circonférentiel en résolvant                         #
c#                            l'équation de Maxwell                                                              #
c#             d2H/dr2+(1/r)dH/dr-(1/r)*H=sigma*mu0*(dM/dt+dH/dt)  (r allant de 0 à 1)                           #
c#                                                                                                               #
c#################################################################################################################
       implicit REAL*8 (A-H,O-Z)
	logical erreur

      parameter Nmax=200,Nhmax=50
	

      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur
      common/para1/Hk,Hz,psi,xnu,xmu0,pi,pi2
	common/eMax/frequence,sigma,xMs,aa,courant,nformeI
      common/paraRK/deltaT,nRK



	dimension ai(Nmax),bi(Nmax),ci(Nmax),di(Nmax)
	dimension Hi(Nmax+2),dUfi(Nmax+2),Hi1(Nmax+2)
      dimension Y(Nmax),Xi(Nmax+2)
              
       
      Rini=0.0d0
	Rfinal=1.0d0
	constante=xmu0*frequence*sigma*aa*aa
      h=(Rfinal-Rini)/(npoint+1)
	x=Rini
	Xi(1)=Rini
	Xi(npoint+2)=Rfinal
	y0=0.0d0
	Hfi0=courant/(pi2*aa*Hk)
	if(nformeI.eq.1) y1=Hfi0*dsin(pi2*t)
	if(nformeI.eq.2) then
	  if(t.lt.5.0d-1) then
	    y1=Hfi0
	  else
	    y1=-Hfi0
	  endif
	endif
	if(nformeI.eq.3) then
	  if(t.lt.5.0d-2) then
	    y1=Hfi0
	  else
	    y1=0.0d0
	  endif
	endif

	do i=1,npoint
	  x=x+h
	  Xi(i+1)=x
	  xx=x*x
	  ai(i)=1.0d0
	  bi(i)=1.0d0/x
        ci(i)=-1.0d0/xx-constante/deltaT
        di(i)=constante*(xMs*dUfi(i+1)/Hk-Hi1(i+1)/deltaT)
c        ci(i)=-1.0d0/xx-constante/deltaT
c        di(i)=constante*(Hi1(i+1)/deltaT)
	enddo

      e1=1.0d0
	e2=0.0d0
	g1=1.0d0
	g2=0.0d0
      ndim=npoint+2
      call resolution(ai,bi,ci,di,Y,e1,e2,g1,g2,y0,y1,ndim)
	do k=1,npoint+2
	  Hi(k)=Y(k)
	enddo

	return
	end

c======================================================================
c=====================================================================
      subroutine MaxwellHz(t,dUz,Hi1,Hi)


c#################################################################################################################
c#                                                                                                               #
c#                 subroutine de calcul du champ magnétique axial en résolvant                                   #
c#                            l'équation de Maxwell                                                              #
c#             d2H/dr2+(1/r)dH/dr=sigma*mu0*(dM/dt+dH/dt)  (r allant de 0 à 1)                           #
c#                                                                                                               #
c#################################################################################################################
       implicit REAL*8 (A-H,O-Z)
	logical erreur

      parameter Nmax=200,Nhmax=50
	

      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur
      common/para1/Hk,Hz,psi,xnu,xmu0,pi,pi2
	common/eMax/frequence,sigma,xMs,aa,courant,nformeI
      common/paraRK/deltaT,nRK
	common/Seb/Hz0,omegaHz,dephas



	dimension ai(Nmax),bi(Nmax),ci(Nmax),di(Nmax)
	dimension Hi(Nmax+2),dUz(Nmax+2),Hi1(Nmax+2)
      dimension Y(Nmax+2),Xi(Nmax+2)
              
       
      Rini=0.0d0
	Rfinal=1.0d0
	constante=xmu0*frequence*sigma*aa*aa
      h=(Rfinal-Rini)/(npoint+1)
	x=Rini
	Xi(1)=Rini
	Xi(npoint+2)=Rfinal
	y0=0.0d0
      y1=Hz0*sin(omegaHz*t+dephas)
	do i=1,npoint
	  x=x+h
	  Xi(i+1)=x
	  ai(i)=1.0d0
	  bi(i)=1.0d0/x
        ci(i)=-constante/deltaT
        di(i)=constante*(xMs*dUz(i+1)/Hk-Hi1(i+1)/deltaT)
	enddo

      e1=0.0d0
	e2=1.0d0
	g1=1.0d0
	g2=0.0d0
	ndim=npoint+2
      call resolution(ai,bi,ci,di,Y,e1,e2,g1,g2,y0,y1,ndim)
	do k=1,npoint+2
	  Hi(k)=Y(k)
	enddo

	return
	end
c========================================================================
c=======================================================================




c======================================================================
c======================================================================

      subroutine resolution(ai,bi,ci,di,Y,e1,e2,g1,g2,y0,y1,ndim)
       implicit REAL*8 (A-H,O-Z)
	parameter Nmax=200 
      logical erreur

c==============================================================c
c                                                              c
c subroutine de résolution de l'équation différentielle        c
c                                                              c
c          (ai*y" + bi*y'+ ci*y = di)                          c
c  avec condition aux limites  e1*y(1)+e2*y'(1)= y0            c
c                              g1*y(ndim)+g2*y'(ndim)= y1      c
c                                                              c
c     ndim= nombre de points = dimension de la matrice         c
c   ai, bi, ci et di sont des tableuax de dimension ndim       c
c                                                              c
c==============================================================c
 
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur

      dimension Y(Nmax+2),ai(Nmax),bi(Nmax),ci(Nmax),di(Nmax),
     &A(Nmax+2,Nmax+2),B(Nmax+2),YY(Nmax+2),pi(Nmax),qi(Nmax),ri(Nmax)  

      h2=h*h
      do k=1,ndim
	  do m=1,ndim
	    A(k,m)=0.0d0
	  enddo
	  B(k)=0.0d0
	enddo


	do j=1,ndim-2
	  pi(j)=2.d0*ai(j)-ci(j)*h2
	  qi(j)=-ai(j)-bi(j)*h/2.0d0
	  ri(j)=-ai(j)+bi(j)*h/2.0d0
	enddo



      do k=2,ndim-1
	 A(k,k)=pi(k-1)
	 A(k,k+1)=qi(k-1)
       A(k,k-1)=ri(k-1)
	 B(k)=-di(k-1)*h2
 	enddo

      B(1)=y0*h
	B(ndim)=y1*h
	A(1,1)=e1*h-e2
	A(1,2)=e2
	A(ndim,ndim-1)=-g2
	A(ndim,ndim)=g1*h+g2

      call resolmat(ndim,A,B,YY)
	do i=1,ndim
	Y(i)=YY(i)
	enddo
	return
	end




c==============================================================================================
c==============================================================================================

      subroutine resolmat(nn,A,B,V)
c==============================================================================================

c subroutine de résolution de l'équation AV=B où A=matrice(nn,nn), B=vecteur(nn) et V=vecteur sollution(nn)
       implicit REAL*8 (A-H,O-Z)
     	 logical erreur

      Parameter Nmax=200
	dimension A(Nmax+2,Nmax+2),B(Nmax+2),V(Nmax+2),nul(Nmax+2)
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur

c  verification des pivots:si un pivot est nul on effectue la permutation de lignes
 	do i=1,nn
	  nul(i)=0

        aii=A(i,i)
	  abaii=dabs(aii)
   	  if(abaii.lt.0.0001) then
	    k=nn
  1       continue
          if(k.eq.i) then
		  nul(i)=1
		  goto 2
	    endif
          aki=A(k,i)
	    abaki=dabs(aki)
	    if(abaki.lt.0.0001) then
		  k=k-1
		  goto 1
		endif
		do n=1,nn 
	      temp=A(i,n)
	      A(i,n)=A(k,n)
	      A(k,n)=temp
	    enddo
	    temp=B(i)
	    B(i)=B(k)
	    B(k)=temp
 2	    continue
	  endif

	enddo



	nn1=nn-1

	do k=1,nn1

        if(nul(i).eq.1) goto 4
c	  akk=A(k,k)
c	  abakk=abs(akk)
c	  if(abakk.lt.0.0001) goto 4

	  do i=k+1,nn
	    fac=A(i,k)/A(k,k)
	    do j=k+1,nn
	      A(i,j)=A(i,j)-fac*A(k,j)
	    enddo
	    B(i)=B(i)-fac*B(k)
	  enddo
 4	  continue
	enddo

      V(nn)=B(nn)/A(nn,nn)
	do k=nn1,1,-1
	  som=0.0
	  do j=k+1,nn
	    som=som+A(k,j)*V(j)
	  enddo
	  V(k)=(B(k)-som)/A(k,k)
	enddo
      aaa=A(nn,nn)

100   format('erreur, collonne',I3,' nulle dans la matrice A dans la  
     & subroutine resolmat')

 3	return
	end
c=========================================================================================================
c=========================================================================================================

  
  
  
  
  
  
  
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
  
  
      subroutine interpole(Xi,Yi,alfa,betta,gamma,etta,nsurf)
 
 
       implicit REAL*8 (A-H,O-Z)
	logical erreur

	parameter Nmax=200
      dimension Xi(Nmax+2),Yi(Nmax+2),alfa(Nmax+1),
     &betta(Nmax+1),gamma(Nmax+1),etta(Nmax+1),A(4,4),B(3,3),sauve(4)
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur

c=================================================================================
c
c   subroutine d'interpolation entre les points de la courbe Yi=f(Xi)
c   
c   entre les points i et i+1  y=alfa(i)*x**3+betta(i)*x**2+gamma(i)*x+etta(i)
c
c les coefficiens alfa(i), ... etta(i) de l'intervalle [i,i+1] sont déterminés par interpolation sur les 4 points i-1, i, i+1 et i+2 
c
c attention: ne pas confondre alpha et alfa
c=================================================================================
      Yi(npoint+3)=Yi(npoint+2)*(1-delR/20.0d0)
	Xi(npoint+3)=Xi(npoint+2)+delR/20.0d0
      npoint=npoint+nsurf
      do n=1,npoint+1
	  alfa(n)=0.0d0
	  betta(n)=0.0d0
	  gamma(n)=0.0d0
	  etta(n)=0.0d0
	enddo
      do i=2,npoint

        netape=0
        x1=Xi(i-1)
	  x2=Xi(i)
	  x3=Xi(i+1)
	  x4=Xi(i+2)
	  y1=Yi(i-1)
	  y2=Yi(i)
	  y3=Yi(i+1)
	  y4=Yi(i+2)
  
construction de la matrice 4x4 de Vandermonde
        xx1=1.0d0
	  xx2=1.0d0
	  xx3=1.0d0
	  xx4=1.0d0
        do n=1,4
	    A(1,n)=xx1
	    A(2,n)=xx2
	    A(3,n)=xx3
	    A(4,n)=xx4
	    xx1=xx1*x1
	    xx2=xx2*x2
	    xx3=xx3*x3
	    xx4=xx4*x4
	  enddo
calcul du déterminant de la matrice A
  1	  do n=1,3
	    do k=1,3
	      B(n,k)=A(n+1,k+1)
	    enddo
	  enddo

        detA=0.0d0
	  xneg=1.0d0

        do k=1,4

          detB=B(1,1)*(B(2,2)*B(3,3)-B(3,2)*B(2,3))
	    detB=detB-B(1,2)*(B(2,1)*B(3,3)-B(3,1)*B(2,3))
	    detB=detB+B(1,3)*(B(2,1)*B(3,2)-B(3,1)*B(2,2))

          detA=detA+xneg*A(1,k)*detB
	    xneg=-xneg
	  
	    if (k.eq.4) goto 2
	    do n=1,3
	      B(n,k)=A(n+1,k)
	    enddo
  2	    continue
	  enddo

	  if(netape.eq.1) goto 3
	  if(netape.eq.2) goto 4
	  if(netape.eq.3) goto 5
	  if(netape.eq.4) goto 6
        det=detA

calcul du coefficient etta
        do m=1,4
	    sauve(m)=A(m,1)
	  enddo
	  A(1,1)=Y1
	  A(2,1)=Y2
	  A(3,1)=Y3
	  A(4,1)=Y4
        

	  netape=1
	  goto 1
 3      continue
        xetta=detA

calcul du coeficient gamma
        do m=1,4
	    A(m,1)=sauve(m)
	    sauve(m)=A(m,2)
	  enddo
	  A(1,2)=Y1
	  A(2,2)=Y2
	  A(3,2)=Y3
	  A(4,2)=Y4

	  netape=2
	  goto 1
 4      continue
        xgamma=detA
calcul du coeficient betta
        do m=1,4
	    A(m,2)=sauve(m)
	    sauve(m)=A(m,3)
	  enddo
	  A(1,3)=Y1
	  A(2,3)=Y2
	  A(3,3)=Y3
	  A(4,3)=Y4

	  netape=3
	  goto 1
 5      continue
        xbetta=detA

calcul du coeficient alfa
        do m=1,4
	    A(m,3)=sauve(m)
	    sauve(m)=A(m,4)
	  enddo
	  A(1,4)=Y1
	  A(2,4)=Y2
	  A(3,4)=Y3
	  A(4,4)=Y4

	  netape=4
	  goto 1
 6      continue
        xalfa=detA

        alfa(i)=xalfa/det
	  betta(i)=xbetta/det
	  gamma(i)=xgamma/det
	  etta(i)=xetta/det


 
	enddo      

      np1=npoint+1
      alfa(1)=alfa(2)
	betta(1)=betta(2)
	gamma(1)=gamma(2)
	etta(1)=etta(2)
      alfa(np1)=alfa(npoint)
	betta(np1)=betta(npoint)
	gamma(np1)=gamma(npoint)
	etta(np1)=etta(npoint)

      Yi(npoint+3)=0.0d0
      npoint=npoint-nsurf
      Xi(npoint+3)=0.0d0
	return
	end




c============================================================================
c
       subroutine integrale(xini,xfin,Xi,alph,betta,gamma,etta,somme)
c
c==============================================================================
c
c   calcul du champ moyen sur une section du cylindre délimitée par les rayons xini et xfin
c la moyenne se fait par integration de (2*pi* rayon * champ) puis divisée par la surface pi* (xfin**2-xini**2)
c le champ est interpolé par les parametres alpha(i), betta(i) et gamma(i)
c            (alpha(i)*x**2+betta(i)*x+gamma(i)) entre chaque deux points Xi(i) et Xi(i+1)
c   l'intervalle d'integration est [xini,xfin] (Xi est le vecteur correspondant aux noeuds de Rini à Rfinal)
c     le résultat est Hmoyen
c     nombre total de noeuds =npoint+2 (point Rini et Rfinal compris)
c==========================================================================================       
       implicit REAL*8 (A-H,O-Z)
	parameter Nmax=200
	logical erreur

	dimension Xi(Nmax+2),alph(Nmax+1),betta(Nmax+1),
     &gamma(Nmax+1),etta(Nmax+1)
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur

      inter1=0
	inter2=npoint+1
      somme=0.0d0
	do i=1,npoint+1
	 if(xini.ge.Xi(i) .and. xini.lt.Xi(i+1))  inter1=i
	 if(xfin.gt.Xi(i) .and. xfin.le.Xi(i+1))  inter2=i
	enddo
	do k=inter1,inter2
	 terme1=alph(k)*(Xi(k+1)**4-Xi(k)**4)/4.0d0
	 terme2=betta(k)*(Xi(k+1)**3-Xi(k)**3)/3.0d0
	 terme3=gamma(k)*(Xi(k+1)**2-Xi(k)**2)/2.0d0
	 terme4=etta(k)*(Xi(k+1)-Xi(k))
	 somme=somme+terme1+terme2+terme3+terme4
	enddo
	terme1=alph(inter1)*(xini**4-Xi(inter1)**4)/4.0d0
	terme2=betta(inter1)*(xini**3-Xi(inter1)**3)/3.0d0
	terme3=gamma(inter1)*(xini**2-Xi(inter1)**2)/2.0d0
	terme4=etta(inter1)*(xini-Xi(inter1))
	surplus1=terme1+terme2+terme3+terme4

	terme1=alph(inter2)*(Xi(inter2+1)**4-xfin**4)/4.0d0
	terme2=betta(inter2)*(Xi(inter2+1)**3-xfin**3)/3.0d0
	terme3=gamma(inter2)*(Xi(inter2+1)**2-xfin**2)/2.0d0
	terme4=etta(inter2)*(Xi(inter2+1)-xfin)
	surplus2=terme1+terme2+terme3+terme4
	somme=somme-surplus1-surplus2
	surf=(xfin**2-xini**2)/2.0d0
	return
	end

      





	subroutine RK4(ti,inst,Ufi,Uz,Ur,Hefi,Hez,Her,dUfi,dUz,dUr)

c  subroutine de calcul des valeurs des composantes de U(r) à l'instant (t + delta t) en partant de leurs valeurs à
c  l'instant (t) avec la méthode de Runge-Kutta 4. les composantes de dU/dt est obtenue par l'équation de LLG.
c  l'argument i correspond au numéro du point sur l'axe radial du fil. 
c
       implicit REAL*8 (A-H,O-Z)
      parameter Nmax=200
      	logical erreur

	common/paraRK/deltaT,nRK
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur
      external eqLLG
	dimension Ufi(Nmax+2),Uz(Nmax+2),Ur(Nmax+2)
	dimension Ufisauve(Nmax+2),Uzsauve(Nmax+2),Ursauve(Nmax+2)
	dimension dUfisauve(Nmax+2),dUzsauve(Nmax+2),dUrsauve(Nmax+2)
	dimension Ufiplus1(Nmax+2),Uzplus1(Nmax+2),Urplus1(Nmax+2)
	dimension UfiA(Nmax+2),UzA(Nmax+2),UrA(Nmax+2)
	dimension UfiB(Nmax+2),UzB(Nmax+2),UrB(Nmax+2)
	dimension UfiC(Nmax+2),UzC(Nmax+2),UrC(Nmax+2)
	dimension dUfi(Nmax+2),dUz(Nmax+2),dUr(Nmax+2)
      dimension Hefi(Nmax+2),Hez(Nmax+2),Her(Nmax+2)

      nRKsauve=nRK
      do j=1,npoint+2
        Ufisauve(j)=Ufi(j)
        Ursauve(j)=Ur(j)
        Uzsauve(j)=Uz(j)
        dUfisauve(j)=dUfi(j)
        dUrsauve(j)=dUr(j)
        dUzsauve(j)=dUz(j)

      enddo
   10 deltaT=1.0d0/inst

      hi=deltaT/nRK

      do n=1,npoint+2
	  Ufiplus1(n)=Ufi(n)
	  Uzplus1(n)=Uz(n)
	  Urplus1(n)=Ur(n)
	enddo

 30	do k=1,nRK

        call eqLLG(Ufi,Uz,Ur,Hefi,Hez,Her,dUfi,dUz,dUr)

        do i=2,npoint+2
	    UfiA(i)=Ufi(i)+dUfi(i)*hi/2.0d0
          UzA(i)=Uz(i)+dUz(i)*hi/2.0d0
	    UrA(i)=Ur(i)+dUr(i)*hi/2.0d0
	    Ufiplus1(i)=Ufiplus1(i)+dUfi(i)*hi/6.0d0
	    Uzplus1(i)=Uzplus1(i)+dUz(i)*hi/6.0d0
	    Urplus1(i)=Urplus1(i)+dUr(i)*hi/6.0d0

	  enddo

	  tihsur2=ti+hi/2.0d0

	  call eqLLG(UfiA,UzA,UrA,Hefi,Hez,Her,dUfi,dUz,dUr)
	  
	  do i=2,npoint+2
	    UfiB(i)=Ufi(i)+dUfi(i)*hi/2.0d0
          UzB(i)=Uz(i)+dUz(i)*hi/2.0d0
	    UrB(i)=Ur(i)+dUr(i)*hi/2.0d0
	    Ufiplus1(i)=Ufiplus1(i)+dUfi(i)*hi/3.0d0
	    Uzplus1(i)=Uzplus1(i)+dUz(i)*hi/3.0d0
	    Urplus1(i)=Urplus1(i)+dUr(i)*hi/3.0d0

	  enddo

	  call eqLLG(UfiB,UzB,UrB,Hefi,Hez,Her,dUfi,dUz,dUr)
	  
	  do i=2,npoint+2
	    UfiC(i)=Ufi(i)+dUfi(i)*hi
          UzC(i)=Uz(i)+dUz(i)*hi
	    UrC(i)=Ur(i)+dUr(i)*hi
	    Ufiplus1(i)=Ufiplus1(i)+dUfi(i)*hi/3.0d0
	    Uzplus1(i)=Uzplus1(i)+dUz(i)*hi/3.0d0
	    Urplus1(i)=Urplus1(i)+dUr(i)*hi/3.0d0

	  enddo
         	
	  tiplush=ti+hi

	  call eqLLG(UfiC,UzC,UrC,Hefi,Hez,Her,dUfi,dUz,dUr)
	  
	  do i=2,npoint+2
	    Ufiplus1(i)=Ufiplus1(i)+dUfi(i)*hi/6.0d0
	    Uzplus1(i)=Uzplus1(i)+dUz(i)*hi/6.0d0
	    Urplus1(i)=Urplus1(i)+dUr(i)*hi/6.0d0

	  enddo

	  do i=2,npoint+2

	    Ufi(i)=Ufiplus1(i)
	    Uz(i)=Uzplus1(i)
	    Ur(i)=Urplus1(i)
        enddo
 	enddo

c   verification de conservation du module de M
      
	do i=1,npoint+2
	  Umodule=dsqrt(Ufi(i)*Ufi(i)+Uz(i)*Uz(i)+Ur(i)*Ur(i))
	  dif=dabs(Umodule-1.0d0)
	  if(dif.gt.1.0d-3) then
	    if(nRK.ge.10000000) then
		  write(6,100) i,ti
	      stop
	    endif
          inst=inst*2
c		nRK=10*nRK
	    do j=1,npoint+2
	      Ufi(j)=Ufisauve(j)
	      Ur(j)=Ursauve(j)
	      Uz(j)=Uzsauve(j)
            dUfi(j)=dUfisauve(j)
            dUr(j)=dUrsauve(j)
            dUz(j)=dUzsauve(j)
		enddo         

         goto 10 
	  endif
	enddo

	do i=1,npoint+2
	  dUfi(i)=(Ufi(i)-Ufisauve(i))/deltaT
	  dUz(i)=(Uz(i)-Uzsauve(i))/deltaT
	  dUr(i)=(Ur(i)-Ursauve(i))/deltaT
	enddo


	nRk=nRKsauve
 100  format('erreur: module de M au point #',I3,' n''est pas conserve a 
     & l''instant t/T=',f6.4)
      return
	end








ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


		subroutine primitive(Ri,alfa,betta,gamma,etta,som)




c  subroutine de calcul de la primitive de la fonction décrite par les parmètres alfa, betta ..etc
       implicit REAL*8 (A-H,O-Z)
	parameter Nmax=200
	logical erreur

      dimension alfa(Nmax+1),betta(Nmax+1),gamma(Nmax+1),etta(Nmax+1)
	dimension som(Nmax+2),Ri(Nmax+2)
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur
      
	do n=1,npoint+2
	som(n)=0.0d0
	enddo
      do n=1,npoint+1
	  Rn=Ri(n)
	  Rn1=Ri(n+1)
	  difsom1=alfa(n)*(Rn1**4-Rn**4)/4.0d0
        difsom2=betta(n)*(Rn1**3-Rn**3)/3.0d0
	  difsom3=gamma(n)*(Rn1**2-Rn**2)/2.0d0
	  difsom4=etta(n)*(Rn1-Rn)
	  difsom=difsom1+difsom2+difsom3+difsom4
	  som(n+1)=som(n)+difsom
	enddo
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


	       subroutine derive(Xi,alfa,betta,gamma,etta,der1,der2)
c
c
c subroutine de calcul des derivés première der1 et deuxième der2 de la fonction définie par:
c    f(xi) = alfa(i)*x**3 + betta(i)*x**2 + gamma(i)*x + etta(i)
c
      implicit REAL*8 (A-H,O-Z)
	parameter Nmax=200
	logical erreur
	dimension Xi(Nmax+2),der1(Nmax+2),der2(Nmax+2)
	dimension alfa(Nmax+1),betta(Nmax+1),gamma(Nmax+1),etta(Nmax+1)
      common/para/const,alpha,delR,h,rapamort,ll,npoint,erreur

      np1=npoint+1
      np2=npoint+2
	do i=1,np1
        x=Xi(i)
	  x2=x*x
        der1(i)=3.0d0*alfa(i)*x2+2.0d0*betta(i)*x+gamma(i)
	  der2(i)=6.0d0*alfa(i)*x+2.0d0*betta(i)	
	enddo

      x=Xi(np2)
	x2=x*x        
      der1(np2)=3.0d0*alfa(np1)*x2+2.0d0*betta(np1)*x+gamma(np1)
	der2(np2)=6.0d0*alfa(np1)*x+2.0d0*betta(np1)	
      return
	end