Télécharger ingase.eso

Retour à la liste

Numérotation des lignes :

  1. C INGASE SOURCE GOUNAND 05/12/21 21:32:14 5281
  2. SUBROUTINE INGASE(MYPGS,IMPR,IRET)
  3. IMPLICIT REAL*8 (A-H,O-Z)
  4. IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM : INGASE
  7. C PROJET : Noyau linéaire NLIN
  8. C DESCRIPTION : Remplit le segment des méthodes d'intégration
  9. C avec des méthodes d'intégration numérique de Gauss
  10. C à une dimension (ordre 1 à 11).
  11. C
  12. C REFERENCES : Numerical recipes (sous-programme gauleg modifié)
  13. C on a recalculé les poids et points de Gauss en REAL*16
  14. C donc avec environ 32 (plutôt 31) chiffres significatifs
  15. C LANGAGE : ESOPE
  16. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  17. C mél : gounand@semt2.smts.cea.fr
  18. C***********************************************************************
  19. C APPELES : INIPG, GCSINO, GCCESY
  20. C APPELE PAR : INPGS
  21. C***********************************************************************
  22. C ENTREES : -
  23. C ENTREES/SORTIES : MYPGS (actif en *MOD)
  24. C SORTIES : -
  25. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  26. C***********************************************************************
  27. C VERSION : v1, 19/10/99, version initiale
  28. C HISTORIQUE : v1, 19/10/99, création
  29. C HISTORIQUE :
  30. C HISTORIQUE :
  31. C***********************************************************************
  32. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  33. C en cas de modification de ce sous-programme afin de faciliter
  34. C la maintenance !
  35. C***********************************************************************
  36. -INC CCOPTIO
  37. CBEGININCLUDE SPOGAU
  38. SEGMENT POGAU
  39. CHARACTER*(LNNPG) NOMPG
  40. CHARACTER*(LNTPG) TYPMPG
  41. CHARACTER*(LNFPG) FORLPG
  42. INTEGER NORDPG
  43. REAL*8 XCOPG(NDLPG,NBPG)
  44. REAL*8 XPOPG(NBPG)
  45. ENDSEGMENT
  46. SEGMENT POGAUS
  47. POINTEUR LISPG(0).POGAU
  48. ENDSEGMENT
  49. CENDINCLUDE SPOGAU
  50. POINTEUR MYPGS.POGAUS
  51. POINTEUR PGCOUR.POGAU
  52. *
  53. INTEGER IMPR,IRET
  54. *
  55. INTEGER DIMSRF
  56. PARAMETER(DIMSRF=1)
  57. REAL*8 XCOR(DIMSRF)
  58. REAL*8 POIDS
  59. *
  60. * Pour les méthodes sur un segment, tous les générateurs sont
  61. * à symétrie centrales
  62. *
  63. * Générateurs pour la méthode de Gauss de degré 1 : GAC1-1-1 :
  64. * - l'origine
  65. REAL*8 X0D1,P0D1
  66. PARAMETER (X0D1=0.D0)
  67. PARAMETER (P0D1=2.D0)
  68. *
  69. * Générateurs pour la méthode de Gauss de degré 3 : GAC1-3-2 :
  70. * - le générateur symétrique (1./\sqrt{3.});(1.)
  71. REAL*8 X1D3,P1D3
  72. PARAMETER (X1D3=0.57735026918962576450914878050195D0)
  73. PARAMETER (P1D3=1.D0)
  74. *
  75. * Générateurs pour la méthode de Gauss de degré 5 : GAC1-5-3 :
  76. * - l'origine
  77. REAL*8 X0D5,P0D5
  78. PARAMETER (X0D5=0.D0)
  79. PARAMETER (P0D5=8.D0/9.D0)
  80. * - le générateur symétrique (\sqrt{3./5.});(5./9.)
  81. REAL*8 X1D5,P1D5
  82. PARAMETER (X1D5=0.77459666924148337703585307995648D0)
  83. PARAMETER (P1D5=5.D0/9.D0)
  84. *
  85. * Générateurs pour la méthode de Gauss de degré 7 : GAC1-7-4 :
  86. * - les générateurs symétriques :
  87. * (\sqrt{\frac{3.-2\sqrt{6./5.}}{7.}});(1./2.+1./(6\sqrt{6./5.}))
  88. * (\sqrt{\frac{3.+2\sqrt{6./5.}}{7.}});(1./2.-1./(6\sqrt{6./5.}))
  89. REAL*8 X1D7,P1D7,X2D7,P2D7
  90. PARAMETER (X1D7=0.33998104358485626480266575910324D0)
  91. PARAMETER (P1D7=0.65214515486254614262693605077800D0)
  92. PARAMETER (X2D7=0.86113631159405257522394648889281D0)
  93. PARAMETER (P2D7=0.34785484513745385737306394922200D0)
  94. *
  95. * Générateurs pour la méthode de Gauss de degré 9 : GAC1-9-5 :
  96. * - l'origine
  97. REAL*8 X0D9,P0D9
  98. PARAMETER (X0D9=0.D0)
  99. PARAMETER (P0D9=128.D0/225.D0)
  100. * - les générateurs symétriques
  101. * ((1./3.)\sqrt{5.-4.\sqrt{5./14.}});(161./450.+13./(180.\sqrt{5./14.}))
  102. * ((1./3.)\sqrt{5.+4.\sqrt{5./14.}});(161./450.-13./(180.\sqrt{5./14.}))
  103. REAL*8 X1D9,P1D9,X2D9,P2D9
  104. PARAMETER (X1D9=0.53846931010568309103631442070021D0)
  105. PARAMETER (P1D9=0.47862867049936646804129151483564D0)
  106. PARAMETER (X2D9=0.90617984593866399279762687829939D0)
  107. PARAMETER (P2D9=0.23692688505618908751426404071992D0)
  108. *
  109. * Générateurs pour la méthode de Gauss de degré 11 : GAC1-11-6 :
  110. * - les générateurs symétriques :
  111. REAL*8 X1D11,P1D11,X2D11,P2D11,X3D11,P3D11
  112. PARAMETER (X1D11=0.23861918608319690863050172168071D0)
  113. PARAMETER (P1D11=0.46791393457269104738987034398956D0)
  114. PARAMETER (X2D11=0.66120938646626451366139959501990D0)
  115. PARAMETER (P2D11=0.36076157304813860756983351383773D0)
  116. PARAMETER (X3D11=0.93246951420315202781230155449399D0)
  117. PARAMETER (P3D11=0.17132449237917034504029614217274D0)
  118. *
  119. INTEGER NOPG
  120. *
  121. * Executable statements
  122. *
  123. IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingase'
  124. *
  125. * Méthode de nom : NCC1-1-2
  126. * Sur un segment : Cubature d'ordre 1 à 2 points
  127. * espace de référence de dimension 1
  128. *
  129. * In INIPG : SEGINI PGCOUR
  130. CALL INIPG('NCC1-1-2','NEWTON-COTES','SEGMENT',
  131. $ 1,2,1,
  132. $ PGCOUR,
  133. $ IMPR,IRET)
  134. IF (IRET.NE.0) GOTO 9999
  135. NOPG=0
  136. XCOR(1)=-1.D0
  137. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET)
  138. XCOR(1)=1.D0
  139. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET)
  140. IF (IRET.NE.0) GOTO 9999
  141. SEGDES PGCOUR
  142. MYPGS.LISPG(**)=PGCOUR
  143. *
  144. * Méthode de nom : NCC1-3-3
  145. * Sur un segment : Cubature d'ordre 3 à 3 points
  146. * espace de référence de dimension 1
  147. *
  148. * In INIPG : SEGINI PGCOUR
  149. CALL INIPG('NCC1-3-3','NEWTON-COTES','SEGMENT',
  150. $ 3,3,1,
  151. $ PGCOUR,
  152. $ IMPR,IRET)
  153. IF (IRET.NE.0) GOTO 9999
  154. NOPG=0
  155. XCOR(1)=0.D0
  156. POIDS=4.D0/3.D0
  157. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET)
  158. XCOR(1)=1.D0
  159. POIDS=1.D0/3.D0
  160. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET)
  161. IF (IRET.NE.0) GOTO 9999
  162. SEGDES PGCOUR
  163. MYPGS.LISPG(**)=PGCOUR
  164. *
  165. * Méthode de nom : GAC1-1-1
  166. * Sur un segment : méthode de Gauss d'ordre 1 à 1 point
  167. * espace de référence de dimension 1
  168. *
  169. * In INIPG : SEGINI PGCOUR
  170. CALL INIPG('GAC1-1-1','GAUSS','SEGMENT',
  171. $ 1,1,1,
  172. $ PGCOUR,
  173. $ IMPR,IRET)
  174. IF (IRET.NE.0) GOTO 9999
  175. NOPG=0
  176. XCOR(1)=X0D1
  177. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D1,IMPR,IRET)
  178. IF (IRET.NE.0) GOTO 9999
  179. SEGDES PGCOUR
  180. MYPGS.LISPG(**)=PGCOUR
  181. *
  182. * Méthode de nom : GAC1-3-2
  183. * Sur un segment : méthode de Gauss d'ordre 3 à 2 points
  184. * espace de référence de dimension 1
  185. *
  186. * In INIPG : SEGINI PGCOUR
  187. CALL INIPG('GAC1-3-2','GAUSS','SEGMENT',
  188. $ 3,2,1,
  189. $ PGCOUR,
  190. $ IMPR,IRET)
  191. IF (IRET.NE.0) GOTO 9999
  192. NOPG=0
  193. XCOR(1)=X1D3
  194. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET)
  195. IF (IRET.NE.0) GOTO 9999
  196. SEGDES PGCOUR
  197. MYPGS.LISPG(**)=PGCOUR
  198. *
  199. * Méthode de nom : GAC1-5-3
  200. * Sur un segment : méthode de Gauss d'ordre 5 à 3 points
  201. * espace de référence de dimension 1
  202. *
  203. * In INIPG : SEGINI PGCOUR
  204. CALL INIPG('GAC1-5-3','GAUSS','SEGMENT',
  205. $ 5,3,1,
  206. $ PGCOUR,
  207. $ IMPR,IRET)
  208. IF (IRET.NE.0) GOTO 9999
  209. NOPG=0
  210. XCOR(1)=X0D5
  211. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D5,IMPR,IRET)
  212. IF (IRET.NE.0) GOTO 9999
  213. XCOR(1)=X1D5
  214. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET)
  215. IF (IRET.NE.0) GOTO 9999
  216. SEGDES PGCOUR
  217. MYPGS.LISPG(**)=PGCOUR
  218. *
  219. * Méthode de nom : GAC1-7-4
  220. * Sur un segment : méthode de Gauss d'ordre 7 à 4 points
  221. * espace de référence de dimension 1
  222. *
  223. * In INIPG : SEGINI PGCOUR
  224. CALL INIPG('GAC1-7-4','GAUSS','SEGMENT',
  225. $ 7,4,1,
  226. $ PGCOUR,
  227. $ IMPR,IRET)
  228. IF (IRET.NE.0) GOTO 9999
  229. NOPG=0
  230. XCOR(1)=X1D7
  231. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D7,IMPR,IRET)
  232. IF (IRET.NE.0) GOTO 9999
  233. XCOR(1)=X2D7
  234. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D7,IMPR,IRET)
  235. IF (IRET.NE.0) GOTO 9999
  236. SEGDES PGCOUR
  237. MYPGS.LISPG(**)=PGCOUR
  238. *
  239. * Méthode de nom : GAC1-9-5
  240. * Sur un segment : méthode de Gauss d'ordre 9 à 5 points
  241. * espace de référence de dimension 1
  242. *
  243. * In INIPG : SEGINI PGCOUR
  244. CALL INIPG('GAC1-9-5','GAUSS','SEGMENT',
  245. $ 9,5,1,
  246. $ PGCOUR,
  247. $ IMPR,IRET)
  248. IF (IRET.NE.0) GOTO 9999
  249. NOPG=0
  250. XCOR(1)=X0D9
  251. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D9,IMPR,IRET)
  252. IF (IRET.NE.0) GOTO 9999
  253. XCOR(1)=X1D9
  254. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D9,IMPR,IRET)
  255. IF (IRET.NE.0) GOTO 9999
  256. XCOR(1)=X2D9
  257. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D9,IMPR,IRET)
  258. IF (IRET.NE.0) GOTO 9999
  259. SEGDES PGCOUR
  260. MYPGS.LISPG(**)=PGCOUR
  261. *
  262. * Méthode de nom : GAC1-11-6
  263. * Sur un segment : méthode de Gauss d'ordre 11 à 6 points
  264. * espace de référence de dimension 1
  265. *
  266. * In INIPG : SEGINI PGCOUR
  267. CALL INIPG('GAC1-11-6','GAUSS','SEGMENT',
  268. $ 11,6,1,
  269. $ PGCOUR,
  270. $ IMPR,IRET)
  271. IF (IRET.NE.0) GOTO 9999
  272. NOPG=0
  273. XCOR(1)=X1D11
  274. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D11,IMPR,IRET)
  275. IF (IRET.NE.0) GOTO 9999
  276. XCOR(1)=X2D11
  277. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D11,IMPR,IRET)
  278. IF (IRET.NE.0) GOTO 9999
  279. XCOR(1)=X3D11
  280. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P3D11,IMPR,IRET)
  281. IF (IRET.NE.0) GOTO 9999
  282. SEGDES PGCOUR
  283. MYPGS.LISPG(**)=PGCOUR
  284. *
  285. * Normal termination
  286. *
  287. IRET=0
  288. RETURN
  289. *
  290. * Format handling
  291. *
  292. *
  293. * Error handling
  294. *
  295. 9999 CONTINUE
  296. IRET=1
  297. WRITE(IOIMP,*) 'An error was detected in subroutine ingase'
  298. RETURN
  299. *
  300. * End of subroutine INGASE
  301. *
  302. END
  303.  
  304.  
  305.  
  306.  

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