c      cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      cc          fitcry routine                                 ccc
c      cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       subroutine fitcry( aa , z , gamma , tem , rhom , qcrst ) 
       implicit real*8 (a-h,o-z)
       real*8 mue,i,j,k,l,ip,jp,kp,lp
       dimension a (0:4),b (4),e(0:4),f(4),i (0:4),j (4),p(0:4),q(4)
     &          ,ap(0:4),bp(4)            ,ip(0:4),jp(4)
       common/wa/ n
c
       mue  = aa/z
       rho  = rhom*mue
       rho6 = rho/1.d6
       t8   = tem/1.d8
       pi   = 3.1415927d0
c
c      -- read parameters ----
       if     ( idint(aa) .eq. 4  ) 
     &    then 
                                   include 'cry-he.par'
       elseif ( idint(aa) .eq. 12 ) 
     &    then 
                                   include 'cry-c.par'
       elseif ( idint(aa) .eq. 16 ) 
     &    then 
                                   include 'cry-o.par'
       elseif ( idint(aa) .eq. 20 ) 
     &    then 
                                   include 'cry-ne.par'
       elseif ( idint(aa) .eq. 24 ) 
     &    then
                                   include 'cry-mg.par'
       elseif ( idint(aa) .eq. 28 ) 
     &    then
                                   include 'cry-si.par'
       elseif ( idint(aa) .eq. 32 ) 
     &    then
                                   include 'cry-s.par'
       elseif ( idint(aa) .eq. 40 ) 
     &    then
                                   include 'cry-ca.par'
       elseif ( idint(aa) .eq. 56 ) 
     &    then
                                   include 'cry-fe.par'
       else
           call errcall( 3 )
           return
       endif
c
       u     = 2.d0*pi*( log10(rho) - 3.d0 )/9.d0
           v =  alpha0 + alpha1 *gamma**(-1.d0/3.d0)
     &                 + alpha2 *gamma**(-2.d0/3.d0) 
     &                 + alpha3 *gamma**(-1.d0)
c
           w =  beta0  + beta1  *gamma**(-1.d0/3.d0)
     &                 + beta2  *gamma**(-2.d0/3.d0) 
     &                 + beta3  *gamma**(-1.d0)
c
c     -------F(u,171)-------
       fb  =  a( 0 )/2.d0  + c*u +d   
      do 100 m = 1, 4
         fm  = float( m)
         fb  =  fb + a( m )*cos( fm*u ) + b( m )*sin( fm*u )
 100  continue
c
c     -------F(u,5000)-------
       ft  =  e( 0 )/2.d0  + g*u +h   
      do 110 m = 1, 4
         fm  = float( m)
         ft  =  ft + e( m )*cos( fm*u ) + f( m )*sin( fm*u )
 110  continue
c
c     -------G(u,171)-------
       gb    =  i( 0 )/2.d0  + k*u +l   
      do 120 m = 1, 4
         fm  = float( m)
         gb  =  gb + i( m )*cos( fm*u ) + j( m )*sin( fm*u )
 120  continue
c
c     -------G(u,5000)-------
       gt    =  p( 0 )/2.d0  + r*u +s   
      do 130 m = 1, 4
         fm  = float( m)
         gt  =  gt + p( m )*cos( fm*u ) + q( m )*sin( fm*u )
 130  continue
c
c
           vp= alphap0 + alphap1*gamma**(-1.d0/3.d0)
     &                 + alphap2*gamma**(-2.d0/3.d0) 
     &                 + alphap3*gamma**(-1.d0)
c
           wp= betap0  + betap1 *gamma**(-1.d0/3.d0)
     &                 + betap2 *gamma**(-2.d0/3.d0) 
     &                 + betap3 *gamma**(-1.d0)
c
c     -------F'(u,171)-------
       fp    =  ap( 0 )/2.d0  + cp*u +dp   
      do 140 m = 1, 4
         fm  = float( m)
         fp  =  fp + ap( m )*cos( fm*u ) + bp( m )*sin( fm*u )
 140  continue
c
c     -------G'(u,171)-------
       gp    =  ip( 0 )/2.d0  + kp*u +lp   
      do 150 m = 1, 4
         fm  = float( m)
         gp  =  gp + ip( m )*cos( fm*u ) + jp( m )*sin( fm*u )
 150  continue
c
c     -------- F_lat(u,Gamma) & G_lat(u,Gamma) -----------
      flat   = ( 1.d0 - v )*fb + v*ft
      glat   = ( 1.d0 - w )*gb + w*gt
c     -------- F_phn(u,Gamma) & G_phn(u,Gamma) -----------
      fphn   = vp*fp
      gphn   = wp*gp
c      
            dn = float( n )
        snthtw = 0.2325d0 
            cv = 0.5d0 + 2.d0*snthtw
            ca = 0.5d0
           cvp = 1.d0  - cv
           cap = 1.d0  - ca
c
        vband = 3.56d-2*z*(rho6/aa)**(1.d0/3.d0)/t8
c       --------------------------------------------------
            if ( vband .ge. 150 ) 
     &        then
                  fband = 0.d0
              else           
                  fband  = exp(-2.d0*vband)
c                  write(*,*) 'vband & fband ',vband,fband
            endif
c       -------------------------------------------------
       qlat    = 0.5738d0*z**2.d0/aa*t8**6.d0*rho*fband
     &          *(0.5d0*(( cv **2.d0 + ca **2.d0 )
     &               + dn*( cvp**2.d0 + cap**2.d0 ))*flat
     &          - 0.5d0*(( cv **2.d0 - ca **2.d0 )
     &               + dn*( cvp**2.d0 - cap**2.d0 ))*glat)
c
c
       qphn    = 0.5738d0*z**2.d0/aa*t8**6.d0*rho
     &          *(0.5d0*(( cv **2.d0 + ca **2.d0 )
     &               + dn*( cvp**2.d0 + cap**2.d0 ))*fphn
     &          - 0.5d0*(( cv **2.d0 - ca **2.d0 )
     &               + dn*( cvp**2.d0 - cap**2.d0 ))*gphn)
c
       qcrst   = qlat + qphn
       return
       end