c Copyright (C) 2008-2014 Vincent Eymet c c KDISTRIBUTION is free software; you can redistribute it and/or modify c it under the terms of the GNU General Public License as published by c the Free Software Foundation; either version 3, or (at your option) c any later version. c KDISTRIBUTION is distributed in the hope that it will be useful, c but WITHOUT ANY WARRANTY; without even the implied warranty of c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the c GNU General Public License for more details. c You should have received a copy of the GNU General Public License c along with KDISTRIBUTION; if not, see c subroutine findroots_3rd(a,b,c,d,nr,x) implicit none include 'max.inc' include 'param.inc' c c Purpose: to find ONLY ONE root of a specified c third-order polynomia: c a*x³+b*x²+c*x+d=0 c c Inputs: c + a,b,c,d: coefficients of the third-degree polynomia c c Outputs: c + nr: number of solutions found c + x: array of solutions c double precision a,b,c,d,x(1:3),z(1:3) integer nr,i double precision p,q,delta double precision u,v,y(1:3) double precision epsilon c label integer strlen character*(Nchar_mx) label label='subroutine findroots_3rd' epsilon=1.0D-5 c Debug c write(*,*) 'a=',a c write(*,*) 'b=',b c write(*,*) 'c=',c c write(*,*) 'd=',d c Debug p=-1.0D+0/3.0D+0*(b/a)**2+c/a q=1.0D+0/27.0D+0*(b/a)*(2.0D+0*(b/a)**2-9.0D+0*c/a)+d/a delta=q**2+4.0D+0/27.0D+0*p**3 c Debug c write(*,*) 'p=',p c write(*,*) 'q=',q c write(*,*) 'delta=',delta c Debug if (delta.gt.0.0D+0) then nr=1 call powalpha((-q+dsqrt(delta))/2.0D+0,1.0D+0/3.0D+0,u) call powalpha((-q-dsqrt(delta))/2.0D+0,1.0D+0/3.0D+0,v) z(1)=u+v c Debug c write(*,*) 'u=',u,u**3,(-q+dsqrt(delta))/2.0D+0 c write(*,*) 'v=',v,v**3,(-q-dsqrt(delta))/2.0D+0 c write(*,*) 'z0=',z0 c Debug else if (delta.eq.0.0D+0) then nr=2 z(1)=3.0D+0*q/p z(2)=-3.0D+0*q/(2.0D+0*p) else if (delta.lt.0.0D+0) then nr=3 do i=1,3 z(i)=2.0D+0*dsqrt(-p/3.0D+0)* & dcos(1.0D+0/3.0D+0* & dacos(-q/2.0D+0*dsqrt(27.0D+0/(-p**3))) & +2.0D+0*i*pi/3.0D+0) enddo endif do i=1,nr x(i)=z(i)-b/(3.0D+0*a) y(i)=a*x(i)**3+b*x(i)**2+c*x(i)+d c if (dabs(y(i)).gt.epsilon) then c write(*,*) 'warning from findroots_3rd :' c write(*,*) 'x(',i,')=',x(i),' a*x³+b*x²+c*x+d=',y(i) c stop c endif enddo c Debug c write(*,*) 'ok findroots_3rd' c Debug return end subroutine dichotomy(nu1,nu2,k1,k2,k0,z1,z2,nu) implicit none include 'max.inc' double precision nu1,nu2,k1,k2,k0 double precision z1,z2 double precision nu double precision S,err double precision x1,x2,x,y1,y2,y double precision epsilon,epsz,errz,minz,signz integer it logical failed c label integer strlen character*(Nchar_mx) label label='subroutine dichotomy' epsz=1.0D-8 minz=z1 if (z2.lt.z2) then minz=z2 endif errz=(dabs(z1)-dabs(z2))/minz signz=z1/z2 c Debug c write(*,*) 'nu1=',nu1,' nu2=',nu2 c write(*,*) 'k1=',k1,' k2=',k2 c write(*,*) 'k0=',k0 c write(*,*) 'z1=',z1,' z2=',z2 c Debug if (nu1.eq.nu2) then write(*,*) 'Error from routine ',label(1:strlen(label)),':' write(*,*) 'nu1=',nu1 write(*,*) 'nu2=',nu2 stop endif if (((k1.lt.k2).and.((k0.lt.k1).or.(k0.gt.k2))).or. & ((k1.ge.k2).and.((k0.gt.k1).or.(k0.lt.k2)))) then write(*,*) 'Error from routine ',label(1:strlen(label)),':' write(*,*) 'k0=',k0,' is out of bounds:' write(*,*) 'k1=',k1 write(*,*) 'k2=',k1 stop endif x1=nu1 y1=k1-k0 x2=nu2 y2=k2-k0 err=dabs(y2) epsilon=dabs(k2-k1)/1.0D+3 if (epsilon.lt.1.0D-10) then epsilon=1.0D-10 endif c Debug c write(*,*) 'epsilon=',epsilon, ' err=',err c Debug if (k0.eq.k1) then nu=nu1 goto 123 endif if (k0.eq.k2) then nu=nu2 goto 123 endif it=0 failed=.false. do while (err.gt.epsilon) it=it+1 x=(x1+x2)/2.0D+0 y=S(z1,z2,nu1,nu2,k1,k2,x)-k0 if (y*y1.gt.0.0D+0) then x1=x y1=y else x2=x y2=y endif err=dabs(y) if (it.gt.Niter_max) then c write(*,*) 'Error from routine ',label(1:strlen(label)),':' c write(*,*) 'Too many iterations' c write(*,*) 'nu1=',nu1 c write(*,*) 'nu2=',nu2 c write(*,*) 'k1=',k1 c write(*,*) 'k2=',k2 c write(*,*) 'k0=',k0 c write(*,*) 'z1=',z1 c write(*,*) 'z2=',z2 c write(*,*) 'x1=',x1,' x2=',x2 c write(*,*) 'y1=',y1,' y2=',y2 c write(*,*) 'epsilon=',epsilon,' err=',err c stop failed=.true. goto 123 endif enddo nu=x if (((nu.lt.nu1).or.(nu.gt.nu2)) c & .and.(errz.le.epsz).and.(signz.lt.0.0D+0) & ) then c linear interpolation ? nu=nu1+(nu2-nu1)/(k2-k1)*(k0-k1) endif 123 continue if (failed) then ! dichotomy procedure fails sometimes (wtf ???) -> linear interpolation too nu=nu1+(nu2-nu1)/(k2-k1)*(k0-k1) endif if ((nu.lt.nu1).or.(nu.gt.nu2)) then write(*,*) 'Error from routine ',label(1:strlen(label)),':' write(*,*) 'nu found=',nu,' is out of bounds:' write(*,*) 'nu1=',nu1 write(*,*) 'nu2=',nu2 write(*,*) 'k1=',k1 write(*,*) 'k2=',k2 write(*,*) 'k0=',k0 write(*,*) 'z1=',z1 write(*,*) 'z2=',z2 stop endif return end subroutine powalpha(x,alpha,y) implicit none double precision x,alpha,y if (x.ge.0.0D+0) then y=x**alpha else y=-(-x)**alpha endif return end