Télécharger moca.eso

Retour à la liste

Numérotation des lignes :

  1. C MOCA SOURCE CHAT 05/01/13 01:46:52 5004
  2. subroutine moca
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8(A-H,O-Z)
  5. -INC CCOPTIO
  6. -INC SMLREEL
  7. pointeur mlree4.mlreel, mlree5.mlreel
  8. character*4 mot(1)
  9. segment itra
  10. integer ipoi(m),npara,ibon(m)
  11. endsegment
  12. segment itrav
  13. real*8 a(n,n),b(n),c(n),bb(n,n)
  14. integer is(n)
  15. endsegment
  16. data MOT /'POID'/
  17. * lecture de la valeur de depart des paramatres
  18. call lirobj ('LISTREEL',mlpar,1,iretou)
  19. mlreel=mlpar
  20. segact mlreel
  21. nparr = prog(/1)
  22. * lecture de la valeur theorique
  23. call lirobj('LISTREEL',mlreel,1,iretou)
  24. if(ierr.ne.0) return
  25. segact mlreel
  26. npmes=prog(/1)
  27. * lecture de la valeur de depart
  28. call lirobj('LISTREEL',mlree1,1,iretou)
  29. if(ierr.ne.0) return
  30. segact mlree1
  31. if(mlree1.prog(/1).ne.npmes) then
  32. write(6,*)' erreur 0 npmes mlree1.prog', npmes,mlree1.prog(/1)
  33. call erreur (19)
  34. return
  35. endif
  36. m=100
  37. segini itra
  38. * write(6,*) ' nb de parametres ',nparr
  39. do 1 i=1,nparr
  40. * lecture des sensibilité par rapport au n parametres
  41. call lirobj('LISTREEL',mlree2,1,iretou)
  42. if(ierr.ne.0) then
  43. write(6,*) ' erreur 1'
  44. call erreur (19)
  45. return
  46. endif
  47. segact mlree2
  48. do iy=1,mlree2.prog(/1)
  49. if( mlree2.prog(iy).ne.0.d0) go to 10
  50. enddo
  51. go to 1
  52. 10 continue
  53. ibon(i)=1
  54. if(mlree2.prog(/1).ne.npmes) then
  55. write(6,*)' erreur 2 npmes mlree2.prog', npmes,mlree2.prog(/1)
  56. call erreur (19)
  57. return
  58. endif
  59. npara=npara+1
  60. ipoi(npara)=mlree2
  61. 1 continue
  62. if( npara.eq.0) then
  63. write(6,*) ' erreur 3 npara ',npara
  64. call erreur (19)
  65. return
  66. endif
  67. * lecture de l'option POIDS
  68. call lirmot(mot,1,itrou,0)
  69. if( itrou.eq.0) then
  70. jg=npmes
  71. isup=1
  72. segini mlree5
  73. do io=1,npmes
  74. mlree5.prog(io)=1.d0
  75. enddo
  76. else
  77. * lecture du poids
  78. call lirobj('LISTREEL',mlree5,1,iretou)
  79. if(ierr.ne.0)return
  80. segact mlree5
  81. isup=0
  82. endif
  83. n = npara
  84. segini itrav
  85. do 3 i=1,npmes
  86. do 4 j=1,npara
  87. mlree3=ipoi(j)
  88. x1=mlree3.prog(i)*2. *mlree5.prog(i)*mlree5.prog(i)
  89. do 5 k=1,npara
  90. mlree4=ipoi(k)
  91. xk=mlree4.prog(i)
  92. a(j,k)=a(j,k)+x1*xk
  93. 5 continue
  94. b(j) = b(j) + x1 * ( mlreel.prog(i)-mlree1.prog(i))
  95. 4 continue
  96. 3 continue
  97. *
  98. ** appel a l'inversion
  99. *
  100. * write(6,*) ' avant inversion'
  101. * write(6,*) ( a(1,it),it=1,npara)
  102. * write(6,*) ( a(2,it),it=1,npara)
  103. * write(6,*) ( a(3,it),it=1,npara)
  104.  
  105. diamax=0.
  106. do io=1,npara
  107. if( abs(a(io,io)).gt.diamax) diamax=abs(a(io,io))
  108. enddo
  109. eps = diamax / 1.e10
  110. call inver (A,npara,icrit,c,is,eps)
  111. if(icrit.ne.0) then
  112. call erreur(21)
  113. return
  114. endif
  115. * write(6,*) ' apres inversion'
  116. * write(6,*) ( a(1,it),it=1,npara)
  117. * write(6,*) ( a(2,it),it=1,npara)
  118. * write(6,*) ( a(3,it),it=1,npara)
  119.  
  120. * write(6,*) ' b '
  121. * write(6,*) ( b(it),it=1,npara)
  122. if(isup.eq.1) then
  123. segsup mlree5
  124. else
  125. segdes mlree5
  126. endif
  127. jg=npara
  128. segini mlree5
  129. do 50 io=1,npara
  130. mlree5.prog(io)=0.d0
  131. do 6 iu=1,npara
  132. mlree5.prog(io)=mlree5.prog(io)+a(io,iu)*b(iu)
  133. 6 continue
  134. 50 continue
  135. * write(6,*) 'valeur trouvées dans moca'
  136. * write(6,*) (mlree5.prog(io),io=1,npara)
  137. segdes mlreel,mlree1
  138. do io=1,npara
  139. mlree3=ipoi(io)
  140. segdes mlree3
  141. enddo
  142. jg=nparr
  143. segini mlreel
  144. ia=1
  145. mlree2=mlpar
  146. do ib=1,nparr
  147. if(ibon(ib).eq.1)then
  148. prog(ib)=mlree5.prog(ia) + mlree2.prog(ib)
  149. ia=ia+1
  150. else
  151. prog(ib)=mlree2.prog(ib)
  152. endif
  153. enddo
  154. segdes mlreel,mlree2
  155. segsup mlree5
  156. call ecrobj('LISTREEL',mlreel)
  157. segsup itrav,itra
  158. return
  159. end
  160.  
  161.  
  162.  
  163.  
  164.  
  165.  
  166.  
  167.  
  168.  

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