Télécharger la plaquette Cast3M


fr en

Télécharger funsho.eso

Retour à la liste

Numérotation des lignes : 

funsho

  1. C FUNSHO    SOURCE    BECC      09/12/07    21:15:35     6579
  2.       subroutine funsho(nordpo, Tmaxcv, acv_r, Rgas_r,
  3.      &     P_r, T_r, u_r, logr,
  4.      &     T, P, rho, u, d)
  5. C
  6. C     INPUT
  7. C
  8. C     nordpo           = order of polynomial for cp and cv (see also
  9. C                        Tmaxcv)
  10. C
  11. C     Tmaxcv           = maximum temperature for cv polynomial expansion
  12. C                        cv(T) = cv(Tmaxcv) if T > Tmaxcv
  13. C
  14. C     acv_r            = to compute cv and ether for r (and rs);
  15. C                        vector such that
  16. C                        cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
  17. C                        ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
  18. C
  19. C     Rgas_r           = gas constant for r and (rs)
  20. C
  21. C     p_r, T_r, u_r    = pressure, temperature and velocities for r
  22. C
  23. C     logr             =  .true. right travelling shock
  24. C                      =  .false. left travelling shock
  25. C
  26. C     T                =  temperature for rs
  27. C
  28. C     OUTPUT
  29. C
  30. C     p, u
  31. C     rho, d           = pressure, velociy, density and
  32. C                        shock speed for rs
  33. C                        NB. for a non-er shock, d = u + c
  34. C                            with c = speed sound in rs
  35. C
  36. c      implicit none
  37.       integer nordpo
  38.       real*8  rindr
  39.       real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
  40.      &     , P_r, T_r, u_r
  41.      &     , T, P, rho, u, d
  42.       real*8 un_r, rho_r
  43.      &     , et, cv, et_r, cv_r, ht, ht_r, b, c, csi, b2
  44.       logical logr
  45. C
  46.       if (logr) then
  47.          rindr = 1.0D0
  48.       else
  49.          rindr = -1.0D0
  50.       endif
  51.       un_r = rindr * u_r
  52. C
  53.       rho_r = P_r / (Rgas_r * T_r)
  54.       call prith1(nordpo, acv_r, Tmaxcv, T_r, et_r, cv_r)
  55.       ht_r = et_r + (Rgas_r * T_r)
  56.       call prith1(nordpo, acv_r, Tmaxcv, T, et, cv)
  57.       ht = et + (Rgas_r * T)
  58. C
  59.       b = ((ht - ht_r) / Rgas_r) + (0.5d0 * (T_r - T))
  60.       c = T_r
  61.       b2 = b * b
  62.       if (((c * T) .lt. (1d-4 * b2)) .and. (b .lt. 0.0D0)) then
  63.          csi = (0.5d0 * c / b2) - (0.125d0 * (T * ((c / b2)**2))) +
  64.      &        (0.0625d0 *  (T*T * ((c / b2)**3)))
  65.          csi = csi * abs(b)
  66.       else
  67.          csi = sqrt(b2 + (c * T)) + b
  68.          csi = csi / T
  69.       endif
  70. C
  71.       rho = rho_R * csi
  72.       P = Rgas_r * rho * T
  73.       u = abs (2.0D0 * (ht - ht_r) *
  74.      &     (csi - 1.0D0) / (csi + 1.0D0))
  75.       u = un_r + (sign(1.0D0,(csi - 1.0D0)) * sqrt(u))
  76. C
  77.       if (T .lt. T_r) then
  78.          d = (cv + Rgas_r) / cv
  79.          d = d * Rgas_r * T
  80.          d = u + sqrt(d)
  81.       else
  82.          if (abs (T - T_r) .lt. (1d-5 * T_r)) then
  83.             d = 2.0D0 * (cv_r + Rgas_r) / cv_r
  84.             d = d * Rgas_r * T_r * csi * csi / (csi + 1.0D0)
  85.          else
  86.             d = 2.0D0 * (ht - ht_r) * csi * csi / ((csi * csi) - 1.0D0)
  87.          endif
  88.          d = abs(d)
  89.          d = un_r + sqrt(d)
  90.       endif
  91. C
  92.       u = rindr * u
  93.       d = rindr * d
  94. C
  95.       return
  96.       end
  97.  
  98.  

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