Télécharger dvboul.eso

Retour à la liste

Numérotation des lignes :

dvboul
  1. C DVBOUL SOURCE KICH 22/09/22 21:15:01 11465
  2. SUBROUTINE DVBOUL(s,centre,ray,nbrobl,nombre,ib,igau)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. c ***********************************************************
  6. c calcul de la plus petite hypersphere, dimension espace < 6
  7. c algorithme inspiré de miniball
  8. c ***********************************************************
  9. integer lisr(nbrobl)
  10. real*8 centre(nbrobl-1),rayon(nbrobl-1)
  11. real*8 s(nbrobl-1,nombre)
  12.  
  13. klisp = 0
  14. 10 continue
  15. lisr(1) = 1
  16. * write(6,*)'advblis1',(lisr(ju),ju = 1,nbrobl)
  17. do iu = 1,nbrobl-1
  18. centre(iu) = s(iu,1)
  19. enddo
  20. ray2 = 0.D0
  21. call dvdmax(s,centre,ray2,lisr,nbrobl,nombre,indx)
  22. * write(6,*)'dvboul1',klisp,((centre(ju)/1.e6),ju=1,nbrobl-1),indx
  23. if(indx.gt.0) then
  24. lisr(2) = lisr(1)
  25. lisr(1) = indx
  26. do iu = 1,nbrobl-1
  27. centre(iu) = (s(iu,lisr(1)) + s(iu,lisr(2)))*0.5D0
  28. rayon(iu) = s(iu,lisr(1)) - centre(iu)
  29. enddo
  30. ray2 = rayon(1)*rayon(1) + rayon(2)*rayon(2) +
  31. &2.D0*rayon(3)*rayon(3)
  32. if (nbrobl.gt.4) then
  33. ray2 = ray2 + 2.D0*(rayon(4)*rayon(4) + rayon(5)*rayon(5))
  34. endif
  35. ray = dsqrt(ray2)
  36. else
  37. * write(6,*) 'boulA',klisp
  38. return
  39. endif
  40. klisp = klisp + 1
  41. * verification
  42. * write(6,*)'advblis2',(lisr(ju),ju = 1,nbrobl)
  43. call dvdmax(s,centre,ray2,lisr,nbrobl,nombre,indx)
  44. * write(6,*)'dvboul2',klisp,((centre(ju)/1.e6),ju=1,nbrobl-1),indx
  45. * nouvel essai 2D
  46. if(indx.gt.0) then
  47. lisr(2) = lisr(1)
  48. lisr(1) = indx
  49. do iu = 1,nbrobl-1
  50. centre(iu) = (s(iu,lisr(1)) + s(iu,lisr(2)))*0.5D0
  51. rayon(iu) = s(iu,lisr(1)) - centre(iu)
  52. enddo
  53. ray2 = rayon(1)*rayon(1) + rayon(2)*rayon(2)
  54. &+2.D0*rayon(3)*rayon(3)
  55. if (nbrobl.gt.4) then
  56. ray2 = ray2 + 2.D0*(rayon(4)*rayon(4) + rayon(5)*rayon(5))
  57. endif
  58. ray = dsqrt(ray2)
  59. klisp = 2
  60. else
  61. * write(6,*) 'boulB',klisp
  62. return
  63. endif
  64. klis0 = 0
  65. NDIM = 2
  66. ray2max = ray2
  67.  
  68. 20 continue
  69. * verification
  70. call dvdmax(s,centre,ray2,lisr,nbrobl,nombre,indx)
  71.  
  72. if(klisp.gt.nombre/2) then
  73. write(6,*) 'pb ', ib,igau
  74. return
  75. endif
  76.  
  77. if (ray2max.gt.0) then
  78. delray = dabs((ray2 - ray2max)/ray2max)
  79. else
  80. delray = 0.d0
  81. endif
  82. if(ib.eq.230.and.igau.eq.6.and..false.) then
  83. write(6,*) 'oul',klisp,klis0,ndim,indx,ray2,delray
  84. write(6,*) ((centre(ju)/1.e6),ju=1,nbrobl-1)
  85. endif
  86. ray = dsqrt(ray2)
  87. if(indx.eq.0.or.delray.lt.1.D-10) then
  88. * write(6,*) 'boulklisp',klisp,ndim
  89. return
  90. else
  91. if (ray2.gt.ray2max) ray2max = ray2
  92. do ll = ndim,2,-1
  93. lisr(ll) = lisr(ll-1)
  94. enddo
  95. lisr(1) = indx
  96. klisp = klisp + 1
  97. if(klisp-1.eq.klis0+ndim) then
  98. ndim = ndim + 1
  99. if (ndim.eq.nbrobl+1) then
  100. c interr(1) = ib
  101. c interr(2) = igau
  102. write(6,*) 'arret boules',klisp, ib,igau,klis0,ndim
  103. call erreur(5)
  104. * call erreur(-366)
  105. return
  106. endif
  107. klis0 = klisp - 1
  108. nd = ndim
  109. else
  110. endif
  111. call dvhyp2(nbrobl,nombre,s,nd,lisr,iok,r2,centre)
  112. if(ib.eq.230.and.igau.eq.6.and..false.) write(6,*) 'R2',r2
  113. ray2 = r2
  114.  
  115. goto 20
  116. endif
  117.  
  118. RETURN
  119. END
  120.  
  121.  
  122.  
  123.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales