home *** CD-ROM | disk | FTP | other *** search
/ Programmer Power Tools / Programmer Power Tools.iso / progjrn / pj_6_1.arc / CHILC.ARC / CHILC.FOR next >
Encoding:
Text File  |  1987-08-31  |  40.3 KB  |  576 lines
  1. C./ ADD NAME=CHILC
  2. C     PROGRAM TO CALCULATE, PRINT, AND PLOT TNE MAGNETIC                
  3. C     SUSCEPTIBILITY OF HSP METALS.                                     
  4. C                                                                       
  5.       IMPLICIT REAL*8 (A-H,P-Z)
  6.       REAL FUN(200,5),NAM(1000),NBOUND(10,2)
  7.       DIMENSION S(1300,3),TM(2,2600),T(2),C(4),S1(3),TE(2,4),TO(2,5)
  8.       COMMON/B/ CRIT, DE ,NTI                                       
  9.       COMMON/A/NA, NS, KL, ISIG                                     
  10.       COMMON/C/EPS, DLE                                             
  11.       COMMON/D/ S12, S11                                            
  12.       REAL*8 KL,KAPPA,KAPRM                                         
  13. C                                                                   
  14.       PI=3.1415927D0                                                
  15.       open (21,FILE='FOR21.DAT')
  16.       open (5,FILE='CHILC.OUT')
  17.  7    READ (21,4100) KODE 
  18.  400  FORMAT(' ENTER RESTART KODE ',I5,/)
  19.       print   400, KODE                                                 
  20.       IF (KODE.LT.1) KODE = 1                                           
  21.       IF (KODE.GT.6) STOP                                               
  22.       GO TO (2, 1, 5, 6, 22, 35), KODE                                  
  23.  2    READ (21,600)  NTI, CRIT, DT, CN, I, NI, NR                       
  24.  500  FORMAT (' ENTER NTI, CRIT, DT, CN, I, NI, NR',I5,3F10.5,3I5,/)    
  25.       print   500, NTI, CRIT, DT, CN, I, NI, NR                         
  26.  600  FORMAT (I5, 3F10.5, 3I5 )                                         
  27. C     NTI IS NUMBER OF SIMPSON INTERVAL IN INTEGRATION FOR TM.          
  28. C     CRIT IS WIDTH (IN KT) OF FERMI FUNCTION USED.                     
  29. C     DT IS RANGE (IN KT) ABOVE FERMI LEVEL WHERE TM CALCULATED.        
  30. C     CN CONTROLS SUBDIVIDING OF SIMPSON INTERVALS IN E INTEGRATION.    
  31. C     IF CN = .0, ENTERED VALUES OF I, NI, NR ARE USEL.                 
  32. C     I AND NI CONTROL RRNGE OF SUBDIVISION, FROM L=I TO L=I+2*NI.      
  33. C     NR IS THE NUMWER OF SUBDIVISIONS PER INTERVAL, SDE=DE/(2NR).      
  34. C     IF NTI = 0, STANDARD VALUES OF ABOVE ARE ENTERED, SEE BELOW.      
  35.       IF (NTI.LT.0) GO TO 7                                             
  36.       IF (NTI.GT.0) GO TO 1                                             
  37.       NTI = 20                                                          
  38.       CRIT = 10.                                                        
  39.       DT = 10.                                                          
  40.       CN = .0                                                           
  41.       I = 0                                                             
  42.       NI = 0                                                            
  43.       NR = 1                                                            
  44. C                                                                       
  45. C     SET UP S VECTORS                                                  
  46.  1    READ (21,2000)  NE, DELTA, EC, BET, ISIG                          
  47.  1000 FORMAT (' ENTER NE, DELTA, EC, BETA, ISIG ',I5,3F10.5,I5,/)       
  48.       print   1000, NE, DELTA, EC, BET, ISIG                            
  49. C     DELTA, EC, AND BET=BETA=B/A ARE BENNETT-FALICOV BAND PARAMETERS.  
  50. C     NE IS THE NUMBER OF SIMPSON INTERVALS, DE = (EMAX-EMIN)/(2NE).    
  51. C     E2/(2DE) AND E3/(2DE) MUST BE INTEGERS, LATER IN THE              
  52. C     PROGRAM, NE BECOMES 2NE+1.                                        
  53. C     IN THIS SECTION IF ISIG GT 1, GET TRACE AND S1,S2,S;              
  54.  2000 FORMAT (I5, 3F10.5, I5)                                           
  55.       IF (NE.LE.0) GO TO 2                                              
  56.       B2 = BET**.66666666D0                                             
  57.       B4 = B2**2                                                        
  58.       RP3 = (B4 + 2./B2)/3.D0                                           
  59.       BE2 = BET**2                                                      
  60.       E3 = EC + DELTA/3.D0                                              
  61.       E2 = EC - DELTA/3.D0                                              
  62.       EMIN = .0                                                         
  63.       IF (E2.LT..0) EMIN = E2                                           
  64.       EMAX = .0                                                         
  65.       IF (E3.GT.0.) EMAX = E3                                           
  66.       EMID = E2                                                         
  67.       IF (EMID.EQ.EMIN) EMID = .0                                       
  68.       IF (EMID.EQ.EMAX) EMID = E3                                       
  69.       EL = .0                                                           
  70.       IF (EC.EQ.0.) GO TO 4                                             
  71.       EL=8.*EC**3/((((4.*RP3/(3.*B4))-1.)*EC**2+DELTA**2/9.D0)*27.*BE2) 
  72.       IF (BET.EQ.1.) GO TO 4                                            
  73.  21   KAPPA = EL*(EL-EC-DELTA/3.)*(EL-EC+DELTA/3.)                      
  74.       KAPRM = 2.*EL*(EL-EC) + ((EL-EC)**2-(DELTA/3.)**2)                
  75.       H3 = RP3*EL - 2.*EC/(3.*B2)                                       
  76.       DEL = -(KAPPA-H3**3)/(KAPRM-3.*RP3*H3**2)                         
  77.       EL = EL + DEL                                                     
  78.       IF (DABS(DEL/(EL-EMID)).GT..001) GO TO 21                         
  79.  4    print  1500, EL, E2, E3                                           
  80.  1500 FORMAT (' EL =',E12.4,' E2 = ',E12.4,' E3 = ',E12.4,/ )           
  81.       DE = (EMAX - EMIN)/(2.*NE)                                        
  82.       NE = 2*NE + 1                                                     
  83.  9    IF (NE.GT.1300) GO TO 90                                          
  84.       IT=0                                                              
  85.       NIT=0                                                             
  86.       NB=0                                                              
  87.       IF (CN.EQ.0.) GO TO 101                                           
  88. C     PICK SUBDIVISION PARAMETERS                                       
  89.       I=0                                                               
  90.       NI=0                                                              
  91.       NR=1                                                              
  92.  101  C(1)=(EMAX-EMID)/DE                                               
  93.       C(2)=(EMID-EMIN)/DE                                               
  94.       C(3)=(EL-EMID)/DE                                                 
  95.       TEMP=DABS(C(3))                                                   
  96.       IF (C(1).GT.C(2)) GOTO 104                                        
  97.       IF (C(1).LT.1.) GOTO 11                                           
  98.       IF (C(1).GE.CN) GOTO 104                                          
  99.       NR = CN/C(1)+1.                                                   
  100.       NI = C(1)+.1                                                      
  101.       I = NE-2*NI                                                       
  102.       IF (I.GE.1) GOTO 107                                              
  103.       NI = NI+(I-1)/2                                                   
  104.       I = 1                                                             
  105.       GOTO 107                                                          
  106.  11   IF (TEMP.LE.2.) GOTO 95                                           
  107.       EMAX=EMAX - 2.*DE                                                 
  108.       NE = NE-2                                                         
  109.       IF (CN.EQ.0.) GOTO 109                                            
  110.       C(1) = C(1)-2.                                                    
  111.       C(3) = C(3)+2.                                                    
  112.       TEMP = TEMP-2.                                                    
  113.       I = NE-2                                                          
  114.       NI=TEMP/4.                                                        
  115.       IF (NI.LT.2) NI=2                                                 
  116.       NR = CN                                                           
  117.       GOTO 107                                                          
  118.  104  IF (C(2).LT.1.) GOTO 12                                           
  119.       IF (C(2).GT.CN) GOTO 107                                          
  120.       NR = CN/C(2)+1.                                                   
  121.       NI = C(2)+.1                                                      
  122.       I = 1                                                             
  123.       GOTO 107                                                          
  124.  12   IF (TEMP.LE.2.) GOTO 95                                           
  125.       EMIN = EMIN+2.*DE                                                 
  126.       NE = NE-2                                                         
  127.       IF (CN.EQ.0.) GOTO 109                                            
  128.       C(2) = C(2)-2.                                                    
  129.       C(3) = C(3)-2.                                                    
  130.       TEMP = TEMP-2.                                                    
  131.       I = 1                                                             
  132.       NI = TEMP/4.                                                      
  133.       IF (NI.LT.2) NI=2                                                 
  134.       NR = CN                                                           
  135.  107  IF (TEMP.GE.CN.OR.CN.EQ.0.) GOTO 109                              
  136.       IF (TEMP.LT..5/CN) GOTO 109                                       
  137.       NRT = CN/TEMP+.999                                                
  138.       NIT = IFIX(SNGL(TEMP+1.))+1                                       
  139.       IT = -1-2*(IFIX(SNGL(TEMP))/4)                                    
  140.       IF (C(3).LT.0.) IT = 2-2*NIT-IT                                   
  141.       IT = IFIX(SNGL(C(2)+.1))+IT                                       
  142.       IF (NRT.GT.NR) NR = NRT                                           
  143.       IF (IT.LT.1) IT = 1                                               
  144.       IF (I.EQ.0) I = IT                                                
  145.       IF (I.LE.IT) GOTO 102                                             
  146.       NI = NI+(I-IT)/2                                                  
  147.       I = IT                                                            
  148.  102  IF (I+2*NI.LT.IT+2*NIT) NI = (IT-I)/2+NIT                         
  149.  109  print  1501, I, NI, NR, NB                                        
  150. C                                                                       
  151. C     ENERGY INTEGRATION                                                
  152.       E = EMIN                                                          
  153.       DO 111 L = 1,NE                                                   
  154.       DO 111 K = 1,3                                                    
  155.  111  S(L,K)=.0                                                         
  156.       IF (ISIG.GT.1) print  1501, I, NI, NR, NB                         
  157.  1501 FORMAT(' I, NI, NR, NB = ',4I4,/)                                 
  158.       NB = 1                                                            
  159.       IF (I.EQ.1)  NB = NR                                              
  160.       IF (ISIG.GT.1) print  1501, I, NI, NR, NB                         
  161.       SDE = DE/NB                                                       
  162.       DLE = .0                                                          
  163.       CALL SS( E, SDE, S1, DELTA, EC, BET)                              
  164.       NEM1= NE-1                                                        
  165.       DO 3 J = 2, NEM1, 2                                               
  166.       NB = 1                                                            
  167.       IF (J.GT.I.AND.J.LT.I+2*NI) NB = NR                               
  168.       IF (ISIG.GT.1) print  1502, I, NI, NR, NB, J                      
  169.  1502 FORMAT (5I4)                                                      
  170.       SDE = DE/NB                                                       
  171.       DO 112 K = 1,3                                                    
  172.  112  S(J-1,K) = S(J-1,K)+S1(K)/(2.*NB)                                 
  173.       IF (ISIG.GT.1) print  6000, E, ( S(J-1,K),K = 1,3)                
  174.       NOSC = 1                                                          
  175.       A = -1.+1./NB                                                     
  176.       NB2 = NB*2                                                        
  177.       DO 115 M = 1,NB2                                                  
  178.       E = E+SDE                                                         
  179. C     IF (J.EQ.NE-1.AND.M.EQ.2*NB) E = EMAX                             
  180. C     IF (DABS(E-EMID).LT..5*SDE.AND.M.EQ.2*NB) E = EMID                
  181. C     DLE = .0                                                          
  182.       EPS = E-EL                                                        
  183.       IF (EPS.GT.-.5*SDE.AND.EPS.LE..5*SDE) DLE = SDE                   
  184.       IF (DLE.EQ.0.D0) GOTO 114                                         
  185.       DO 113 K = 1,3                                                    
  186.  113  TM(1,K) = S1(K)                                                   
  187.       IF (DLE.NE.0.D0.AND.DABS(E-EMID).LT..5*SDE) GOTO 17               
  188.  114  CALL SS( E, SDE, S1, DELTA, EC, BET)                              
  189.       IF (DLE.EQ.0.D0) GOTO 15                                          
  190. C                                                                       
  191. C     CALCULATE CONTRIBUTIONNOF SINGULARITY AT E = EL                   
  192.       DO 13 K = 1,3                                                     
  193.  13   C(K) = S1(K)                                                      
  194.       DLE = .0D0                                                        
  195.       CALL SS( E+SDE, SDE, S1, DELTA, EC, BET)                          
  196.       IF (DABS(E-EMID).GT..1*SDE) GOTO 119                              
  197.       S1(1) = S1(1)+S11                                                 
  198.       S1(2) = S1(2)+S12                                                 
  199.  119  TEMP = EPS/SDE                                                    
  200.       DO 14 K = 1,3                                                     
  201.       IF (MOD(M,2).NE.0) S1(K)=(2.*S1(K)+2.*TM(1,K)+(3.D0*TEMP          
  202.      1 *DLOG((1.+TEMP)/(1.-TEMP))-6.D0)*C(K))/4.D0                      
  203.       IF (MOD(M,2).EQ.0) S1(K) = (S1(K)+TM(1,K)+(3.D0*TEMP              
  204.      1 *DLOG((1.+TEMP)/(1.-TEMP))+.5D0*DLOG((4.-TEMP**2)/(1.-TEMP**2)   
  205.      2 )-6.D0)*C(K))/2.                                                 
  206.  14   CONTINUE                                                          
  207.       GOTO 15                                                           
  208. C                                                                       
  209. C     CALCULATE CONTRIBUTION OF SINGULARITY AT E = EL = EMID            
  210.  17   DLE = .0D0                                                        
  211.       DTE = DABS(EMID-2.*EC/(3.*RP3*B2))                                
  212.       IF (DTE.GT..5*SDE) DTE = .5*SDE                                   
  213.       IF (DTE.LT..05*SDE) DTE = .05*SDE                                 
  214.       TW3 = 2.**.33333333D0                                             
  215.       TEMP = 3.D0/(2.D0*TW3)-1.                                         
  216.       TEMP3 = (SDE/DTE)**.33333333D0-2.+DTE/SDE+(1.-DTE/SDE)/TW3        
  217.       TEMP1 = 3.D0*TEMP/TEMP3                                           
  218.       TEMP2 = (1.+3.D0*TEMP*(DTE/SDE-2.)/TEMP3)                         
  219.       TEMP3 = -.5D0+3.*TEMP*(1.-DTE/SDE)/TEMP3                          
  220.       CALL SS ( E-2.*SDE, SDE, C, DELTA, EC, BET)                       
  221.       DO 117 K = 1,3                                                    
  222.  117  TM(1,K) = TEMP2*S1(K)+TEMP3*C(K)                                  
  223.       CALL SS ( E+SDE, SDE, S1, DELTA, EC, BET)                         
  224.       CALL SS ( E+2.*SDE, SDE, C, DELTA, EC, BET)                       
  225.       DO 121 K = 1,3                                                    
  226.  121  TM(1,K) = TM(1,K)+TEMP2*S1(K)+TEMP3*C(K)                          
  227.       CALL SS ( E-DTE, DTE, S1, DELTA, EC, BET)                         
  228.       CALL SS ( E+DTE, DTE, C, DELTA, EC, BET)                          
  229.       DO 118 K = 1,3                                                    
  230.  118  S1(K) = TM(1,K)+TEMP1*(S1(K)+C(K))                                
  231. C                                                                       
  232.  15   TEMP = (3.D0+NOSC)/(4.D0*NB)                                      
  233.       IF ( M.EQ.2*NB) TEMP = TEMP/2.                                    
  234.       IF (ISIG.GT.1) print  1503, M, NOSC, A, TEMP                      
  235.  1503 FORMAT (' M, NOSC, A, TEMP ', 2I4, 2E12.4)                        
  236.       DO 116 K = 1,3                                                    
  237.       S(J-1,K) = S(J-1,K)-TEMP*A*(1.-A)*S1(K)                           
  238.       S(J,K) = S(J,K)+TEMP*(1.-A**2)*S1(K)                              
  239.  116  S(J+1,K) = S(J+1,K)+TEMP*A*(1.+A)*S1(K)                           
  240.       IF (ISIG.GT.1) print  6000, E, ( S(J-1,K), K = 1,3)               
  241.       IF (ISIG.GT.1) print  6000, E, ( S(J,K), K = 1,3)                 
  242.       IF (ISIG.GT.1) print  6000, E, ( S(J+1,K), K = 1,3)               
  243.       A = A+1./NB                                                       
  244.   115 NOSC = -NOSC                                                      
  245.       IF (DABS(E-EMID).GT..5*SDE.OR.DABS(EPS).LT..5*SDE) GOTO 3         
  246. C                                                                       
  247. C     STEP FUNCTION CORRECTOR                                           
  248.       TEMP = S11*DSIGN(1.D0,EL-EMID)                                    
  249.       S(J+1,1) = S(J+1,1)-TEMP/(2.*NB)                                  
  250.       S1(1) = S1(1)+TEMP                                                
  251.       TEMP = S12*DSIGN(1.D0,EL-EMID)                                    
  252.       S(J+1,2) = S(J+1,2)-TEMP/(2.*NB)                                  
  253.       S1(2) = S1(2)+TEMP                                                
  254. C                                                                       
  255.  3    CONTINUE                                                          
  256.       DO 120 K = 1,3                                                    
  257.       S(1,K) = 2.*S(1,K)                                                
  258.  120  S(NE,K) = 2.*S(NE,K)                                              
  259. C                                                                       
  260. C     PLOT OR PRINT S; ISIG = 0, PLOT S(M); ISIG > 0. PRINT;            
  261. C     ISIG < 0. PRINT AND PLOT.                                         
  262.       READ (21,4100)  ISIG,M                                            
  263.  3000 FORMAT (' ENTER ISIG,M ',2I5,/)                                   
  264.       print   3000, ISIG,M                                              
  265.  4100 FORMAT ( 2I5 )                                                    
  266.       IF (ISIG) 5,5,6
  267.  5    IF (M.LT.1.OR.M.GT.3) GOTO 20                                     
  268.       SMIN = S(1,M)                                                     
  269.       SMAX = SMIN                                                       
  270.       DO 16 L = 2,NE                                                    
  271.       SMIN = DMIN1(SMIN,S(L,M))                                         
  272.  16   SMAX = DMAX1(SMAX,S(L,M))                                         
  273.       print  4200, SMIN, SMAX                                           
  274.  4200 FORMAT (' SMIN =',E12.4,' SMAX =',E12.4/' ENTER SMIN,SMAX '/)     
  275. C     READ (21,4000), SMIN, SMAX                                        
  276.  4000 FORMAT (2F10.5)                                                   
  277. C     CALL PLOT (EMIN,.0,12)                                            
  278. C     CALL PLOT (.0,SMIN,13)                                            
  279. C     CALL PLOT (.0,SMIN,12)                                            
  280. C     CALL PLOT (.0,SMAX,12)                                            
  281. C     CALL TTY                                                          
  282.  6    IF (ISIG.EQ.3) print  5000                                        
  283.  5000 FORMAT(' ',5X,'E',8X,'S(1)',8X,'S(2)',8X,'N(E)')                  
  284.       E = EMIN                                                          
  285.       DO 10 L = 1,NE                                                    
  286.       IF (ISIG.EQ.3) print  6000, E, ( S(L,K), K = 1,3)                 
  287.  6000 FORMAT (F12.6, 3E12.4)                                            
  288. C     IF (ISIG.LE.0.AND.L.EQ.1) CALL PLOT(EMIN,S(L,M),13)               
  289. C     IF (ISIG.LE.0) CALL PLOT ( E, S(L,M), 12)                         
  290. C     IF (ISIG.LT.0) CALL TTY                                           
  291.  10   E = E+DE                                                          
  292. C     CALL TTY                                                          
  293.       IF (ISIG.LE.0) GOTO 5                                             
  294. C                                                                       
  295. C     PREPARE S WITH SIMPSON FACTORS                                    
  296.  20   OSC = DE/3.D0                                                     
  297.       DO 8 L = 1,NE                                                     
  298.       OSC = -OSC                                                        
  299.       SIMP = DE+OSC                                                     
  300.       IF (L.EQ.1.OR.L.EQ.NE) SIMP = DE/3.D0                             
  301.       DO 8 K = 1,3                                                      
  302.  8    S(L,K) = S(L,K)*SIMP                                              
  303. C                                                                       
  304.       LGR = 1                                                           
  305. C     SET UP TM'S                                                       
  306.  22   READ (21,8000)  EFMIN, EFMAX, TRY, A, B                           
  307.  7000 FORMAT (' ENTER  EFMIN, EFMAX, TRY, ABAR, B ',5F10.6,/)           
  308.       print   7000, EFMIN, EFMAX, TRY, A, B                             
  309.  8000 FORMAT (5F10.6)                                                   
  310. C     EFMIN, EFMAX ARE MINIMUM AND MAXIMUM FERMI LEVELS                 
  311. C     TRY IS THE TEMPERATURE (IN ENERGY UNITS)                          
  312. C     ABAR AND B ARE BAND PARAMETERS (NOT THE B/A =BETA)                
  313.       IF (B.EQ.0.) GOTO 1                                               
  314. C     FMUL = 1.71732D-6*A**2/DSQRT(B)                                   
  315.       FMUL = 1.71732D-6*A**2/DSQRT(B)   /7.14                           
  316.       FNE = 1.8769D4*RP3/A**4                                           
  317.       NM = -(EFMIN-EMIN)/DE                                             
  318.       EFMIN = EMIN-NM*DE                                                
  319.       NM = (EFMAX-EMIN)/DE                                              
  320.       EFMAX = EMIN+NM*DE                                                
  321.       NF = (EFMAX-EFMIN)/DE+1.1                                         
  322.       NT = (EFMAX-EMIN+DT*TRY)/DE+2.1                                   
  323.       IF (NT.GT.NE+NF) NT = NE+NF                                       
  324.       IF (NT.GT.2600) GOTO 93                                           
  325.       IF (NT.LT.1) GOTO 35                                              
  326.       EMU = EMIN-EFMAX                                                  
  327.       DO 30 L = 1,NT                                                    
  328.       CALL TN( 1.D0, TRY, EMU, T)                                       
  329.       IF (ISIG.EQ.4) print  6000, EMU, ( T(K), K = 1,2)                 
  330.       DO 25 K = 1,2                                                     
  331.  25   TM(K,L) = T(K)*FMUL                                               
  332.  30   EMU = EMU+DE                                                      
  333.       IF (TRY.NE.0.) GOTO 35                                            
  334.       CALL TONE(2.D0,TO)                                                  
  335.       SM = DSQRT(DE)*FMUL                                               
  336.       DO 32 K = 1,2                                                     
  337.       DO 31 N = 1,4                                                     
  338.  31   TO(K,N) = TO(K,N)*SM                                              
  339.       IF (ISIG.GT.2) print  6000, SM, (TO(K,L), L=1,4)                  
  340.  32   SM = SM/DE                                                        
  341.       OSC = 1.D0/3.D0                                                   
  342.       DO 34 N = 1,4                                                     
  343.       DO 33 K = 1,2                                                     
  344.       TE(K,N) = TO(K,N)/(1.-OSC)                                        
  345.  33   TO(K,N) = TO(K,N)/(1.+OSC)                                        
  346.  34   OSC = -OSC                                                        
  347.       DO 28 K = 1,2                                                     
  348.       TO(K,5) = TO(K,4)                                                 
  349.       TO(K,4) = TO(K,3)                                                 
  350.       TO(K,3) = TO(K,2)                                                 
  351.       TO(K,2) = TO(K,1)+3.D0*TM(K,NT-3)/8.D0                            
  352.       TO(K,1) = 5.D0*TM(K,NT-4)/4.D0                                    
  353.  28   IF (ISIG.GT.2) print  6000, OSC, (TO(K,L), L = 1,5)               
  354.       DO 29 K = 1,2                                                     
  355.       TM(K,NT-3) = TE(K,1)+TM(K,NT-3)/2.                                
  356.       TM(K,NT-2) = TE(K,2)                                              
  357.       TM(K,NT-1) = TE(K,3)                                              
  358.  29   TM(K,NT) = TE(K,4)                                                
  359. C                                                                       
  360. C     PLOT OR PRINT C'S AND CHI = C(3), N(E) = C(4)                     
  361. C     IF ISIG = 0, PLOT C(M)                                            
  362. C     IF ISIG > 0, PRINT ALL C'S                                        
  363. C     IF ISIG < 0, PRINT ALL C'S AND PLOT C(M)                          
  364.  35   READ (21,8800)  ISIG, M,NPR                                       
  365.  8500 FORMAT(' ENTER ISIG, M, NPR',3I5,/)                               
  366.       print   8500, ISIG, M, NPR                                        
  367.  8800 FORMAT (3I5)                                                      
  368.       IF (M.LE.0.OR.M.GT.4) GOTO 22                                     
  369.       IF (ISIG) 36,36,37                                                
  370.  36   print  9000                                                       
  371.  9000 FORMAT (' ENTER CMIN, CMAX '/)                                    
  372.       READ (21,4000)  CMIN, CMAX                                        
  373.       DO 52 K=1,5                                                       
  374.       NBOUND(K,1) = CMIN                                                
  375.  52   NBOUND(K,2) = CMAX                                                
  376.  9700 FORMAT(9X,'EF',4X,'C(1)',7X,'C(2)',7X,'CHI',8X,'N(E)',/)          
  377.  37   IF (ISIG.GE.1) print  9700                                        
  378. C                                                                       
  379. C     COMRUTE C'S                                                       
  380.       EF = EFMIN                                                        
  381.       DO 50 N = 1,NF                                                    
  382.       DO 38 K = 1,4                                                     
  383.  38   C(K) = .0D0                                                       
  384. C                                                                       
  385. C     END CORRECTIONS                                                   
  386.       NTI = 2*NTI                                                       
  387.       CRIT = 2.*CRIT                                                    
  388.       DEN = .0D0                                                        
  389.       IF (DABS(E2).LT.DE) DEN = DE                                      
  390.       CALL SS ( EMIN, DEN, S1, DELTA, EC, BET)                          
  391.       IF (DABS(E2).LT.DE) S12 = S1(2)                                   
  392.       CALL TN(1.D0, TRY, EMIN-EF, T)                                    
  393.       C(1) = -FMUL*T(1)*S12                                             
  394.       DEN = .0D0                                                        
  395.       IF (DABS(E3).LT.DE) DEN = DE                                      
  396.       CALL SS ( EMAX, DEN, S1, DELTA, EC, BET)                          
  397.       IF (DABS(E3).LT.DE) S12 = S1(2)                                   
  398.       CALL TN( 1.D0, TRY, EMAX-EF, T)                                   
  399.       C(1) = C(1)+FMUL*T(1)*S12                                         
  400.       IF (DABS(E2).GT.DE.AND.DABS(E3).GT.DE) GOTO 45                    
  401. C                                                                       
  402. C     END CORRECTIONS IF E2 OR E3 = 0                                   
  403.       CALL TN( 1.D0, TRY, -EF, T)                                       
  404.       C(2) = C(2)-PI*FMUL*T(2)/(3.*RP3*B2)                              
  405.       C(1) = C(1)+FMUL*T(1)*PI*(2.*B4/3.D0+.5D0*(RP3+1./B2))/(RP3*EC)   
  406.       TEMP = -2.D0*FMUL*PI*B4*DE/(3.D0*RP3*EC*DABS(EC))                 
  407.       C(1) = C(1)+TEMP*T(1)                                             
  408.       CALL TN (1.D0, TRY, DSIGN(DE,EC)-EF, T)                           
  409.       C(1) = C(1)+4.D0*TEMP*T(1)                                        
  410.       CALL TN( 1.D0, TRY, 2.*DSIGN(DE,EC)-EF, T)                        
  411.       C(1) = C(1)+TEMP*T(1)                                             
  412.  45   NTI = NTI/2                                                       
  413.       CRIT = CRIT/2.D0                                                  
  414.       NL = NT+N-NF                                                      
  415.       IF (NL.LT.3) GOTO 41                                              
  416.       IF (NL.GT.NE) NL = NE                                             
  417.       DO 40 L = 1,NL                                                    
  418.       NTL = NF-N+L                                                      
  419.       IF (NTL.GT.NT) GOTO 41                                            
  420.       TMM = TM(2,NTL)                                                   
  421.       IF(TRY.NE.0) GOTO 39                                              
  422.       IF (MOD(NL,2).NE.0.AND.NTL.GE.NT-4) TMM = TO(2,NTL-NT+5)          
  423.  39   C(4) = C(4)+S(L,3)*TMM                                            
  424.       DO 40 K = 1,2                                                     
  425.       TMM = TM(K,NTL)                                                   
  426.       IF (TRY.NE.0) GOTO 40                                             
  427.       IF (MOD(NL,2).NE.0.AND.NTL.GE.NT-4) TMM = TO(K,NTL-NT+5)          
  428.  40   C(K) = C(K)+S(L,K)*TMM                                            
  429.  41   DO 42 K = 1,2                                                     
  430.  42   C(3) = C(3)+C(K)                                                  
  431.       C(4) = FNE*C(4)                                                   
  432.       NN=N/NPR
  433.       IF (ISIG.GE.1.AND.NN*NPR.EQ.N) print  9800, EF, ( C(K), K = 1,4)  
  434.  9800 FORMAT (F12.8,4D11.3)                                             
  435.       IF (NN*NPR.EQ.N) NAM((LGR-1)*(NF/NPR)+NN) = C(M)                  
  436.       IF (LGR.NE.1.OR.N.NE.1) GO TO 48                                  
  437.       CMIN = C(M)                                                       
  438.       CMAX = C(M)                                                       
  439.  48   CMIN = DMIN1(CMIN, C(M))                                          
  440.       CMAX = DMAX1(CMAX, C(M))                                          
  441.  49   CONTINUE                                                          
  442.  50   EF = EF+DE                                                        
  443. C     DO 52 K=1,LGR                                                     
  444. C     NBOUND(K,1) = CMIN                                                
  445. C52   NBOUND(K,2) = CMAX                                                
  446.       LGR=LGR+1                                                         
  447.       IF (ISIG.GE.0) GO TO 22                                           
  448.       LGR=LGR-1                                                         
  449.       CALL GRAF (EFMIN,EFMAX,NAM,NBOUND,NF/NPR,LGR  ,FUN)               
  450.       print  9900, CMIN, CMAX                                           
  451.  9900 FORMAT (' CMIN = ', E12.4,' CMAX = ', E12.4)                      
  452. C                                                                       
  453.       IF (ISIG.EQ.2) GO TO 1                                            
  454.       IF (ISIG.EQ.-1) STOP                                              
  455. C                                                                       
  456. C     ERROR MESSAGES                                                    
  457.    90 print  1100                                                       
  458.  1100 FORMAT (' NE TOO LARGE')                                          
  459.       GOTO 1                                                            
  460.  93   print  1200                                                       
  461.  1200 FORMAT(' NT TOO LARGE')                                           
  462.       GOTO 22                                                           
  463.  95   print  1300                                                       
  464.  1300 FORMAT (' NE TOO SMALL')                                          
  465.       GOTO 1                                                            
  466.       END                                                               
  467. C                                                                       
  468.       SUBROUTINE TN(B, TRY, EMU, T )                                    
  469. C                                                                       
  470.       IMPLICIT  REAL *8 (A-H,O-Z)                                       
  471. C     SUBROUTINE TO CALCULATE THE INTEGRAL OF THE FERMI FUNCTION        
  472. C     OF E = B*Z**2, T(1), AND ITS DERIVATIVE WITH RESPECT TO           
  473. C     FERMI LEVEL, T(2), TRY IS THE TEMPERATURE (IN ENERGY UNITS)       
  474. C     AND EMU = E - MU,WHERE MU IS THE FERMI LEVEL,                     
  475.       DIMENSION T(2)                                                    
  476.       COMMON/B/     CRIT, DE ,NT                                        
  477. C                                                                       
  478.       PI = 3.1415927D0                                                  
  479.       IF (DABS(EMU).LT..1*DE) EMU = .0D0                                
  480.       T(1) = .0D0                                                       
  481.       T(2) = .0D0                                                       
  482.       N = NT                                                            
  483.       IF (TRY.EQ.0.) GO TO 40                                           
  484.       BETA = 1.D0/TRY                                                   
  485.       BEMU = BETA*EMU                                                   
  486.       IF (BEMU.GT.CRIT) GO TO 70                                        
  487.       IF (BEMU.LT.-CRIT) GO TO 80                                       
  488.       Z0 = 0.D0                                                         
  489.       ARF = ((-CRIT)/BETA-EMU)/B                                        
  490. C     IF (ARF.GT.0D0)  Z0 = DSQRT(ARF)                                  
  491.       ZL = DSQRT((((CRIT+5.D0)/BETA) -EMU)/B)                           
  492. C                                                                       
  493.       DZ = (ZL-Z0) /(2*N )                                              
  494.       OSC = DZ/3.D0                                                     
  495.       Z = Z0                                                            
  496.       N2P1  = 2*N  + 1                                                  
  497.       DO 50 L = 1, N2P1                                                 
  498.       OSC = -OSC                                                        
  499.       EX = .0                                                           
  500.       SIMP = DZ + OSC                                                   
  501.       IF (L.EQ.1.OR.L.EQ.2*N +1) SIMP = DZ/3.D0                         
  502.       BZ2 = B*Z*Z                                                       
  503.       ARG = BETA * ( EMU + BZ2 )                                        
  504.       IF (DABS(ARG).GT.4.*CRIT) GO TO 50                                
  505.       EX = DEXP(ARG)                                                    
  506.       D = 1.D0 + EX                                                     
  507.       DFDM = BETA*EX/D**2                                               
  508.       T(1) = T(1) + 2.D0*SIMP*BZ2*DFDM                                  
  509.       T(2) = T(2) + SIMP*DFDM                                           
  510.  50   Z = Z + DZ                                                        
  511.       IF (EX.EQ.0D0) RETURN                                             
  512.       TEMP = .5/(EX*BETA*B*(Z-DZ))                                      
  513.       DO 60 L = 1, 2                                                    
  514.       T(L) = T(L) + TEMP                                                
  515.  60   TEMP = BETA*TEMP                                                  
  516.       RETURN                                                            
  517. C                                                                       
  518.  40   IF (EMU.GE.0.D0) RETURN                                           
  519.       T(1) = DSQRT(-EMU/B)                                              
  520.       T(2) = -.5D0*T(1)/EMU                                             
  521.       RETURN                                                            
  522. C                                                                       
  523.  70   IF (BEMU.GT.50) RETURN                                            
  524.       EX = DEXP(-BEMU)                                                  
  525.       TEMP = .5D0*EX*DSQRT(PI/(BETA*B))                                 
  526.       DO 75 L = 1, 2                                                    
  527.       T(L) = TEMP                                                       
  528.  75   TEMP = TEMP*BETA                                                  
  529.       RETURN                                                            
  530. C                                                                       
  531. C                                                                       
  532.  80   T(1) =DSQRT(-EMU/B)                                               
  533.       T(2) = -.5D0*T(1)/EMU                                             
  534.       TEMP = (.5D0*T(2)/EMU)*((PI/BETA)**2)/6.D0                        
  535.       T(1) = T(1) + TEMP                                                
  536.       T(2) = T(2) + 1.5*TEMP/EMU                                        
  537.       RETURN                                                            
  538.       END                                                               
  539. C                                                                       
  540. C                                                                       
  541.       SUBROUTINE TONE(X,T)                                              
  542.       IMPLICIT  REAL *8 (A-H,O-Z)                                       
  543.       DIMENSION T(2,4)                                                  
  544.       X2 = X*X                                                          
  545.       R = DSQRT(X)                                                      
  546.       T(1,1) = X2*R*(X2/9.D0-.2D0 )/3.D0                                
  547.       T(2,1) = X*R*(X2/6.D0-X2/7.D0+1.D0/9.D0-1.D0/6.D0)                
  548.       T(1,2) = X2*R*(-X2/9.D0+X/7.D0+.4D0 )                             
  549.       T(2,2) = X*R*(3.D0*X2/7.D0-2.D0*X/5.D0-2.D0/3.D0-(X+1.)*(X-2.)/2.)
  550.       T(1,3) = X*R*(X*X2/9.D0-2.D0*X2/7.D0-X/5.D0+2.D0/3.D0)            
  551.       T(2,3) = -X*R*(3.D0*X2/7.D0-4.D0*X/5.D0-1.D0/3.D0)                
  552.       T(1,4) = -X2*R*(X2/9.D0-3.D0*X/7.D0+2.D0/5.D0)/3.D0               
  553.       T(2,4) = X*R*(3.D0*X2/7.D0-6.D0*X/5.D0+2.D0/3.D0)/3.D0            
  554.       RETURN                                                            
  555.       END                                                               
  556.       SUBROUTINE GRAF (EFMIN,EFMAX,NAM,NBOUND,NF,LGR,FUN)
  557.       LOGICAL*1 FIELD,EQSMBL,SMBLS(5)
  558.       LOGICAL LEQ,LALLX
  559.       INTEGER GR,DIGY,GRINTX,FSTNX
  560.       REAL *8 LBX,STX,LBY,SCY,EFMIN,EFMAX
  561.       REAL NAM(310),NBOUND(10,2),FUN(NF,LGR)
  562.       DATA FIELD/1H /,EQSMBL/1H*/,SMBLS/1H*,1H+,1H=,1H^,1H-/
  563.       LALLX=.FALSE.
  564.       LEQ=.FALSE.
  565.       DO 22 L=1,LGR
  566.       DO 22 N=1,NF
  567.  22   FUN(N,L)=NAM((L-1)*NF+N)
  568.       LBX=EFMIN
  569.       STX=(EFMAX-EFMIN)/(NF-1)
  570.       LBY=NBOUND(1,1)
  571.       SCY=NBOUND(1,2)
  572.       CALL MAZIL2(FIELD,EQSMBL,LEQ,LGR,NF,FUN,SMBLS,LBX,STX,
  573.      *            LBY,SCY,3,3,10,1,LALLX,0)
  574.       RETURN
  575.       END
  576.