U
    _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   s4   | dkrdS d}t d| d D ]}||9 }q|S d S )Nr      r      )range)nvali r	   S/home/ubuntu/Home/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 r   r	   r   )r   r	   r
   r   %   s     )sumr   )startstopr   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 )%
KineticGas      ?r   c                    sB  | 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| _t|D ]p\}}tdd t|t|D ||d	 D ]@\}	}
| t|t|
| j||	f< | j||	f | j|	|f< qqd
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   d   r   r   c                 S   s   g | ]}|qS r	   r	   r   r   r	   r	   r
   r)   S   s     N      ?   )*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_fracsr@   Zcomplistpq_ranger   pr   qdelta_0bdrQ   rR   rS   r	   r'   r
   __init__2   s4    

".":<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 )Nr,   r   r   r!   )r   rJ   rK   rL   r?   r   r0   rM   rG   rD   rI   rN   rO   r8   r9   r1   rT   r?   rZ   r[   r\   rQ   rR   rS   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 )Nr,   r   r   g      )r   rJ   rK   rL   r?   r   r0   rM   rG   rD   rI   rN   rO   r8   r9   r1   r6   r7   r^   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 )Nr,   r   r   r-   )r   rJ   rK   rL   r?   r   r0   rM   rG   rD   rI   rN   rO   r8   r6   r9   r7   r/   )	rT   r?   rZ   r[   r\   rQ   rR   rS   alpha_T0r	   r	   r
   ra   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&   r5   r   )locZilocr%   dfr	   r
   r)      s     z5KineticGas.get_hard_sphere_radius.<locals>.<listcomp>r"   r!   r   )axis)pd
read_excelospathdirname__file__r   r/   r.   r   meshgridvstack)rT   rU   Zsigma_iZsigma_ijr	   rc   r
   r3      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   )      r!   (z, z$) are non-valid arguments for omega.)r5   r   rJ   pirK   rL   r?   r0   r   r>   r1   r6   r7   
ValueErrorstr)rT   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   rX   rY   r   r	   r
   inner   s    <T
*zKineticGas.A.<locals>.innerr   r   minrT   rX   rY   r   r   r}   r	   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 rw   )r   r6   r7   r   )r   r   r   k)FGr   rX   rY   r   rT   r	   r
   r}      s8    8

X
.Bz!KineticGas.A_prime.<locals>.innerc                    s,   t dt d  |  || fdS )Nr   r   r   r~   )r   r   )r}   rX   rY   r   r	   r
   sum_w   s    z!KineticGas.A_prime.<locals>.sum_wc                    s   t  d t | | dS )Nr   r   r~   r{   )r   r   r	   r
   sum_k   s    z!KineticGas.A_prime.<locals>.sum_kr   r6   r7   r   r   )rT   rX   rY   r   r   r   r	   )	r   r   r}   r   rX   rY   r   rT   r   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 rw   rz   r{   r|   r	   r
   r}      s    8
T
*z'KineticGas.A_tripleprime.<locals>.innerr!   r   r~   r   r	   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 )	Nrp   c                    s     | |d||  S Nro   )r   rv   r   r   rX   rY   rT   r	   r
   r}      s    zKineticGas.H_ij.<locals>.innerc                    s   t |  d |   | dS Nr   r   r   r   r}   rX   rY   r	   r
   sum_r   s    zKineticGas.H_ij.<locals>.sum_rrx   r!   r   r   )rT   rX   rY   ru   r6   r7   r   r   r	   r}   rX   rY   rT   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 )Nrp   c                    s     | |d||  S r   )r   rv   r   r   r	   r
   r}      s    zKineticGas.H_i.<locals>.innerc                    s   t |  d |   | dS r   r   r   r   r	   r
   r      s    zKineticGas.H_i.<locals>.sum_rrx   r   )r7   r6   r   r   )rT   rX   rY   ru   r   r   r	   r   r
   H_i   s    zKineticGas.H_ic                    s<    fddfdd}dt dtd | S )Nc                    s    | | ||  S r   )r   rv   r   )r   rX   rY   rT   r	   r
   r}      s    z"KineticGas.H_simple.<locals>.innerc                    s   t |  d |   | dS r   r   r   r   r	   r
   r      s    z"KineticGas.H_simple.<locals>.sum_rrx   r   r   r~   )rT   rX   rY   r   r   r	   )r   r}   rX   rY   rT   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   r!   ro   rp   r   r   )r6   r8   r9   r   r7   r   r   )rT   rX   rY   r	   r	   r
   rH      s"    $($("2zKineticGas.ac                 C   s   t | j| | j | gS r   )r   r/   rP   )rT   r?   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 }td
D ]}|	|
| ||< qd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< dD ]6}|	j|
| |d d ||< tj||  d!d" q.td
D ]}td#|d  t||  || | _t| j| jg| _d| j| j  | _t|D ]"\}}td$| td|g| _t| j| _| j\| _| _| jt| j | _| j\| _| _t|D ]\}}t|d| g| _| j\| _ | _!t| |d d}t"d| d d| d ft#| _$t|D ]r\}}t%d%d t|t&|D ||d& D ]@\}}| 't(|t(|| j$||f< | j$||f | j$||f< qȐqd'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< q0qt%||D ].\}}tj0||t1|d||t2| d( q|d	krZt||  t3|| d || d  t4dd tj5d)t6t1| j| j d d*d+ qnt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   r!   T   ry   r-   rx   r   )r   r-   )figsizer   )ncolsnrowsfigurewspaceZwidth_ratiosc                 S   s   g | ]}d qS r   r	   r   r	   r	   r
   r)      s     z(KineticGas.plot_test.<locals>.<listcomp>)皙?g)\()r   g333333ÿ)r   gɿblack)colorsalphagffffff?g?r   )r   r   r   r   c                 S   s   g | ]}d qS r   r	   r   r	   r	   r
   r)   .  s     rn   )ZshareyF)ZvisiblezMaking plotzm =c                 S   s   g | ]}|qS r	   r	   r+   r	   r	   r
   r)   L  s     Nr,   )labelcolorz$\frac{\sigma_2}{\sigma_1} = $   )fontsizez$\frac{m_2}{m_1}$gffffff?)titlebbox_to_anchoriX  )dpi):cmget_cmapr   linspacerM   r=   r/   pltr   gsZGridSpecr   Zadd_subplotscaZhlinesrA   setpZget_yticklabelsprintr<   r5   r>   rE   r0   r   r1   r:   r;   r2   r6   r7   rV   r8   r9   rB   rC   rD   rF   rG   rH   rI   rJ   rK   rL   r?   rN   rO   rP   plotroundmaxylimZxlimr   rt   legend	get_titlesavefig)rT   r@   comparesavecmapx_listZs_listZ
sigma_listZm_listfigZgridZaxsr   Zlim_listZs_linemr   xrW   indrX   r   rY   rZ   r[   r\   rQ   rR   rS   sr   r	   r	   r
   	plot_test  s~    





". :<(
(zKineticGas.plot_testN)TF)__name__
__module____qualname__r]   r_   r`   ra   r3   rv   r   r   r   r   r   r   rH   r   r   r	   r	   r	   r
   r    0   s   /	
	
	r    )N)__doc__numpyr   Zscipy.linalglinalgrN   scipy.constants	constantsrK   matplotlib.pyplotpyplotr   matplotlib.cmr   matplotlib.gridspecgridspecr   pandasrf   rh   r   r   r   r   r    r	   r	   r	   r
   <module>   s   	
