Esope source files

   1 : C WEIPA1    SOURCE    CHAT      05/01/13    04:11:29     5004           
   2 : c***********************************************************************        
   3 : c Metodo della  Massima  Verosimiglianza . Soluzione iterativa  della  *        
   4 : c equazione seguente con il metodo delle bisezioni :                   *        
   5 : c                                                                      *        
   6 : c m=N/(N*SUM(sig(i)**m*log(sig(i)))/SUM(sig(i)**m)-SUM(log(sig(i)))) *
   7 : c                                                                      *        
   8 : c***********************************************************************
   9 : c
  10 :       SUBROUTINE WEIPA1(n,sig,volut,rappca,rml,sigzer)
  11 : c
  12 :       IMPLICIT INTEGER(I-N)
  13 :       IMPLICIT REAL*8(A-H,O-Z)
  14 :       dimension sig(n)
  15 : c
  16 : c precisione di calcolo della soluzione                                         
  17 :       eps=1.0e-6                                                                
  18 : c incremento di 'm' per la ricerca dell'intervallo contenente la soluz.         
  19 :       rincr=0.1e+1                                                              
  20 :       rn=float(n)                                                               
  21 : c valore di partenza di 'm'                                                     
  22 :       rml=1.0e-6                                                                
  23 : c valore di partenza della differenza tra i membri dell'equazione               
  24 :       diffa=0.                                                                  
  25 :       do 11 j=1,1000                                                            
  25 :         call emme(n,rml,sig,sum1,diff)                                            
  26 :         diffb=diff                                                                
  27 :         if(diffa*diffb.lt.0.)then                                                 
  28 :           a1=rml1                                                                   
  29 :           b2=rml                                                                    
  30 :           c2=log((b2-a1)/eps)/log(2.)+1.
  31 :           l=idint(c2)                                                                
  32 :           do 3 k=1,l                                                                
  33 :             rml=(b2+a1)/2.                    
  33 :             call emme(n,rml,sig,sum1,diff)           
  34 :             if(diff*diffa.lt.0.)then              
  35 :               b2=rml                
  36 :             else                                       
  37 :               a1=rml                                            
  38 :               diffa=diff                                      
  39 :             endif                                    
  40 :    3      continue                                                                  
  41 :           go to 999                                                                 
  42 :         else                                                                      
  43 :           rml1=rml                                                                  
  44 :           diffa=diffb                                                               
  45 :           rml=rml+rincr                                                             
  46 :         endif                                                                     
  47 :   11  continue
  48 : C                                                                  
  49 :  999  continue
  50 : C                                                                  
  51 :       b1=rn/sum1
  52 : C     
  53 : c calcolo del volume effettivo
  54 :       rkv=(rappca*rml+1.)/(2.*(rml+1.)**2)
  55 :       volef=volut*rkv
  56 : C 
  57 : C sigma zero    
  58 :       sigzer = (volef / b1 ) ** (1./rml)
  59 : C 
  60 : C funzioni gamma
  61 :       uno  = 1.+ 1./rml                                                             
  62 :       rdue = 1.+ 2./rml                                                            
  63 :       gam1 = dgamma(uno)
  64 :       gam2 = dgamma(rdue)
  65 : C     
  66 : C valore medio                                                
  67 :       sigmed=b1**(-1./rml)*gam1                                                 
  68 : C
  69 : C varianza e deviazione standard                                      
  70 :       var   = b1**(-2./rml)*(gam2-gam1*gam1)                                        
  71 :       devst = sqrt(var)                                                           
  72 : C
  73 : *      write(6,13)volut
  74 : *  13  format(/,1x,'VOLUME UTILE                 :',e15.8)
  75 : *      write(6,991)volef
  76 : *  991 format(1x,'VOLUME EFFETTIVO DEL PROVINO :',e15.8)
  77 :       write(6,400)
  78 :  400  format(1x,'Mean stress    ',1x,'  Variance     ')          
  79 :       write(6,500)sigmed,var                                    
  80 :  500  format(2(1x,e15.8),/)                                                       
  81 : c        
  82 :       return
  83 :       end                                                                       
  84 :  
  85 :  
  86 :