B
    _E9                 @   s   d Z ddlZddlmZ ddlmZ ddl	m
Z ddlmZ ddlmZ ddlZddlZdd ZdddZdd Zd	d
 ZG dd dZdS )u  
Author: Vegard G. Jervell
Date: December 2020
Purpose: Implementation of the Chapman-Enskog solutions to the Boltzmann equations for a binary system
        as proposed by Tompson, Tipton and Loyalka, in the paper
        "Chapman–Enskog solutions to arbitrary order in Sonine polynomials III:
        Diffusion, thermal diffusion, and thermal conductivity in a binary, rigid-sphere, gas mixture"
        doi : https://doi.org/10.1016/j.euromechflu.2008.12.002

Requires : numpy, scipy, matplotlib, pandas
    Nc             C   s8   | dkrdS d}xt d| d D ]}||9 }q W |S d S )N)r      r      )range)nvali r   V/Users/vegardjervell/Documents/7_semester/irrev_prosjekt_renskriv/models/kineticgas.pyfac   s    r
   c                sN    d k	r*t  fddt| |d D S t fddt| |d D S d S )Nc             3   s   | ]}| V  qd S )Nr   ).0r   )argsfuncr   r	   	<genexpr>#   s    zsummation.<locals>.<genexpr>r   c             3   s   | ]} |V  qd S )Nr   )r   r   )r   r   r	   r   %   s    )sumr   )startstopr   r   r   )r   r   r	   	summation!   s    "r   c             C   s   | |krdS dS d S )Nr   r   r   )r   jr   r   r	   delta'   s    r   c             C   s0   ddd| d  dd|      t j|d  S )Ng      ?r   r   )npmath	factorial)lrr   r   r	   w-   s    r   c               @   s   e Zd ZddgdfddZdd Zdd Zd	d
 Zdd Zdd Zdd Z	dd Z
dd Zdd Zdd Zdd Zdd Zdd Zd$d!d"Zd#S )%
KineticGasg      ?r   c                sJ  | d}t fdd|D | _t| j| _| j| j | _| || _t	| j| _
| j\| _| _|\| _| _| j\| _| _| j
\| _| _| jd | _d| _|| _t| |d d}td| d d| d ft| _xt|D ]t\}}xjtdd t|t|D ||d	 D ]@\}	}
| t|t|
| j||	f< | j||	f | j|	|f< qW qW d
t t!j"| j t| j  }t#d| d }||tt|d d < t$%| j|}|tt|d d d tt|d d d  \}}}dd|   | j| t | j | j| t | j   | _&|||  | _'| _(| _)d	S )a  
        :param comps (str): Comma-separated list of components, following Thermopack-convention
        :param eos (thermopack): An initialized equation of state, only used to get mole weights
        :param mole_fracs (array-like): list of mole fractions
        :param N (int > 0): Order of approximation.
                            Be aware that N > 10 can be detrimental to runtime. This should be implemented in C++ or Fortran.
        ,c                s   g | ]}   |qS r   )compmoleweightgetcompindex)r   comp)eosr   r	   
<listcomp>=   s    z'KineticGas.__init__.<locals>.<listcomp>)r   r   d   r   r   c             S   s   g | ]}|qS r   r   )r   r   r   r   r	   r"   S   s    Ng      ?   )*splitr   arraymole_weightsr   m0Mget_hard_sphere_radiusZsigmaijdiagsigmaM1M2x1x2m1m2sigma1sigma2sigma12TNarangeemptyfloatA_matr	enumeratezipr   lenaintsqrtconst	Boltzmannzeroslinsolvesoretd_1d0d1)selfcompsr!   
mole_fracsr7   Zcomplistpq_ranger   pr   qdelta_0bdrH   rI   rJ   r   )r!   r	   __init__2   s4    

"0"":<zKineticGas.__init__c             C   s   dt tj| j t | j  }t t| j	}||t
t|d d < t| j	|}|t
t|d d d t
t|d d d  \}}}d| j | j t dtj | | j  | S )Ng      ?r   r   g      ?)r   rA   rB   rC   r6   r   r'   rD   r>   r;   r@   rE   rF   r/   r0   r(   )rK   r6   rQ   rR   rS   rH   rI   rJ   r   r   r	   interdiffusiona   s    ":zKineticGas.interdiffusionc             C   s   dt tj| j t | j  }t t| j	}||t
t|d d < t| j	|}|t
t|d d d t
t|d d d  \}}}d| j | j t dtj | | j  | j| t | j | j| t | j   S )Ng      ?r   r   g      )r   rA   rB   rC   r6   r   r'   rD   r>   r;   r@   rE   rF   r/   r0   r(   r-   r.   )rK   r6   rQ   rR   rS   rH   rI   rJ   r   r   r	   thermal_diffusionj   s    ":(zKineticGas.thermal_diffusionc       	      C   s   dt tj| j t | j  }t t| j	}||t
t|d d < t| j	|}|t
t|d d d t
t|d d d  \}}}dd|   | j| t | j | j| t | j   }t | |gS )Ng      ?r   r   r$   )r   rA   rB   rC   r6   r   r'   rD   r>   r;   r@   rE   rF   r/   r-   r0   r.   r&   )	rK   r6   rQ   rR   rS   rH   rI   rJ   alpha_T0r   r   r	   rW   t   s    "::zKineticGas.alpha_T0c                sX   t tjtd  t fdd|dD }dtj	t
|t|dd }|S )a  
        Get hard-sphere diameters, assumed to be equal to Mie-potential sigma parameter. Gets Mie-parameters from the file mie.xlsx.

        :param comps: Comma seperated list of components
        :return: N x N matrix of hard sphere diameters, where sigma_ij = 0.5 * (sigma_i + sigma_j),
                such that the diagonal is the radius of each component, and off-diagonals are the average diameter of
                component i and j.
        z	/mie.xlsxc                s(   g | ] } j  d  |k d jd qS )r    r,   r   )locZiloc)r   r    )dfr   r	   r"      s    z5KineticGas.get_hard_sphere_radius.<locals>.<listcomp>r   g      ?r   )axis)pd
read_excelospathdirname__file__r   r&   r%   r   meshgridvstack)rK   rL   Zsigma_iZsigma_ijr   )rY   r	   r*      s    	 z!KineticGas.get_hard_sphere_radiusc             C   s   |dkrJ| j |d  d ttjtj | j | j|d    t|| S |dkrd| j	d  tdtj tj | j | j
| j | j   t|| S tdt| d t| d t| d d S )	N)r   r   r   r   )      g      ?(z, z$) are non-valid arguments for omega.)r,   r   rA   pirB   rC   r6   r'   r   r5   r(   r-   r.   
ValueErrorstr)rK   ijr   r   r   r   r	   omega   s
    BFzKineticGas.omegac                s8    fdd}t  d t d  |S )Nc                s  d|  t  d|    d   d|    t d  t d d |    d  t |  t |   t   t | d    t |   t  d |    t d d  t  d |   d d    | d    d |     |     S )N   r   r   r      )r
   )r   )r   rO   rP   r   r   r	   inner   s    ZzKineticGas.A.<locals>.innerr   )r   min)rK   rO   rP   r   r   rm   r   )r   rO   rP   r   r	   A   s    zKineticGas.Ac                s   j d jd  dj  j   j j j  fddfddfdd}td t d  |S )	Nr   c                s  |\}}d| t  d|  |   d|   t d  t d d | |    dd    ||   |   t | |  t | |   t |  t  d |  |   t d d  t  d | |   d d   t | t ||  t |   dd|  d  j|  j | |     d j d |  |   t| j|  t|d    S )Nrk   r   r   r   rl   )r
   r-   r.   r   )r   r   r   k)FGr   rO   rP   r   rK   r   r	   rm      s     Cz!KineticGas.A_prime.<locals>.innerc                s,   t dt d  |  || fdS )Nr   r   )r   )r   rn   )rp   r   )rm   rO   rP   r   r   r	   sum_w   s    z!KineticGas.A_prime.<locals>.sum_wc                s   t  d t | | dS )Nr   )r   )r   rn   )r   )r   rs   r   r	   sum_k   s    z!KineticGas.A_prime.<locals>.sum_kr   )r-   r.   r   rn   )rK   rO   rP   r   r   rt   r   )	rq   rr   rm   r   rO   rP   r   rK   rs   r	   A_prime   s    $
zKineticGas.A_primec          	      sX    d dkrdS  fdd}d d  t  d t d  | S )Nr   r   c                s  d|  t  d|    d d|    t d  t d d |    dd   t |  t |   t   t | d    t |   t  d |    t d d  t  d |   d d    | d    d |     |     S )Nrk   r   r   r   rl   )r
   )r   )r   rO   rP   r   r   r	   rm      s    Zz'KineticGas.A_tripleprime.<locals>.innerg      ?r   )r   rn   )rK   rO   rP   r   r   rm   r   )r   rO   rP   r   r	   A_tripleprime   s    zKineticGas.A_tripleprimec                sv   j j }}|dkr || }}fdd  fdd}d|d   |d   tdtd | }|S )	Nrd   c                s     | |d||  S )Nrc   )ro   rj   )r   r   )rO   rP   rK   r   r	   rm      s    zKineticGas.H_ij.<locals>.innerc                s   t |  d |   | dS )Nr   )r   )r   )r   )rm   rO   rP   r   r	   sum_r   s    zKineticGas.H_ij.<locals>.sum_rrk   g      ?r   )r-   r.   r   rn   )rK   rO   rP   ri   r-   r.   rw   r   r   )rm   rO   rP   rK   r	   H_ij   s    
2zKineticGas.H_ijc                sr   |dkrj j __ fdd  fdd}dtdtd | }|dkrnj j __ |S )Nrd   c                s     | |d||  S )Nrc   )ru   rj   )r   r   )rO   rP   rK   r   r	   rm      s    zKineticGas.H_i.<locals>.innerc                s   t |  d |   | dS )Nr   )r   )r   )r   )rm   rO   rP   r   r	   rw      s    zKineticGas.H_i.<locals>.sum_rrk   r   )r.   r-   r   rn   )rK   rO   rP   ri   rw   r   r   )rm   rO   rP   rK   r	   H_i   s    zKineticGas.H_ic                s<    fddfdd}dt dtd | S )Nc                s    | | ||  S )N)rv   rj   )r   r   )r   rO   rP   rK   r   r	   rm      s    z"KineticGas.H_simple.<locals>.innerc                s   t |  d |   | dS )Nr   )r   )r   )r   )rm   rO   rP   r   r	   rw      s    z"KineticGas.H_simple.<locals>.sum_rrk   r   r   )r   rn   )rK   rO   rP   r   rw   r   )r   rm   rO   rP   rK   r	   H_simple   s    zKineticGas.H_simplec             C   s  |dks|dkr|dkr<| j d | j | j | ||d S |dk rl| jd  | j | j | | |d S |dkr| j d | j | j | ||d S |dk r| jd  | j | j | || d S | j | j | j | ||d S n|dkr0|dkr0| jd | ||d | j| j | ||d  S |dkr`|dk r`| j| j | || d S |dk r|dkr| j| j | | |d S | jd | | | d | j| j | | | d  S d S )Nr   g      ?rc   rd   r   r   )r-   r/   r0   ry   r.   rz   rx   )rK   rO   rP   r   r   r	   r?      s"    $($("2zKineticGas.ac             C   s   t | j| | j | gS )N)r   r&   rG   )rK   r6   r   r   r	   get_alpha_T0  s    zKineticGas.get_alpha_T0TFc          
   C   s  t d}tddd}td}dddg}d| _|d	krtddd
ddddg}tjdd}	t	j
dd|	dddddgd}
dd td
D }x"td
D ]}|	|
| ||< qW dddg}t|d  tjdddddd ntddd}tjdd}	t	j
d
d|	dd}
dd td
D }|	|
d |d< x>dD ]6}|	j|
| |d d ||< tj||  d!d" q4W xBtd
D ]4}td#|d  t||  || | _t| j| jg| _d| j| j  | _x>t|D ]0\}}td$| td|g| _t| j| _| j\| _| _| jt| j | _| j\| _| _xt|D ]\}}t|d| g| _| j\| _ | _!t| |d d}t"d| d d| d ft#| _$xt|D ]v\}}xjt%d%d t|t&|D ||d& D ]@\}}| 't(|t(|| j$||f< | j$||f | j$||f< qW qW d't)t*j+| j, | j  }td| d }||t(t&|d d < t-.| j$|}|t(t&|d d d t(t&|d d d  \}}}dd|   | j | t)| j | j!| t)| j   | _/| j/|||f< qDW qW x<t%||D ].\}}tj0||t1|d||t2| d( qW |d	kr~t||  t3|| d || d  t4dd tj5d)t6t1| j| j d d*d+ qzW tj7d,dd-gd.}tj|8 d*d+ |rtj9|d/d0 d&S )1aR  
        Plot soret coefficient from kinetic gas theory for N'th order approximation with different m2/m1 and sigma2/sigma1 ratios
        :param N (int): Order of approximation
        :param compare (bool): Use same limits as Tipton, Tompson and Loyalka
        :param save (str): Filename to save figure as (defaults to False)
        viridisgMbP?g+?2   )
   r}   r   r   g      ?T   rl   r$   rk   r~   )r~   r$   )figsizer   )ncolsnrowsfigurewspaceZwidth_ratiosc             S   s   g | ]}d qS )Nr   )r   r   r   r   r	   r"      s    z(KineticGas.plot_test.<locals>.<listcomp>)g?g)\()g?g333333ÿ)r   gɿblack)Zcolorsalphagffffff?g?g?)r   r   r   r   c             S   s   g | ]}d qS )Nr   )r   r   r   r   r	   r"   .  s    )r   r   )ZshareyF)ZvisiblezMaking plotzm =c             S   s   g | ]}|qS r   r   )r   r   r   r   r	   r"   L  s    Ng      ?)labelcolorz$\frac{\sigma_2}{\sigma_1} = $   )fontsizez$\frac{m_2}{m_1}$gffffff?)titlebbox_to_anchoriX  )dpi):cmget_cmapr   linspacerD   r4   r&   pltr   gsZGridSpecr   Zadd_subplotscaZhlinesr8   setpZget_yticklabelsprintr3   r,   r5   r<   r'   r   r(   r1   r2   r)   r-   r.   rM   r/   r0   r9   r:   r;   r=   r>   r?   r@   rA   rB   rC   r6   rE   rF   rG   plotroundmaxylimZxlimr   rh   legend	get_titlesavefig)rK   r7   comparesavecmapx_listZs_listZ
sigma_listZm_listfigZgridZaxsr   Zlim_listZs_linemrp   xrN   indrO   r   rP   rQ   rR   rS   rH   rI   rJ   sr   r   r   r	   	plot_test  s~    






"0$:<*
*zKineticGas.plot_testN)TF)__name__
__module____qualname__rT   rU   rV   rW   r*   rj   ro   ru   rv   rx   ry   rz   r?   r{   r   r   r   r   r	   r   0   s   /	
	
	r   )N)__doc__numpyr   Zscipy.linalglinalgrE   scipy.constants	constantsrB   matplotlib.pyplotpyplotr   matplotlib.cmr   matplotlib.gridspecgridspecr   pandasr[   r]   r
   r   r   r   r   r   r   r   r	   <module>   s   	
