介绍

python语言的科研绘图合集，数据来源Hydrogen-diffusion-and-water-rock-reaction

导入python包

import pandas as pd
import glob
import proplot as pplt
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import patheffects
#import ternary
import scipy as scp
import matplotlib as mpl
import matplotlib.colors as mcolors
import matplotlib.cm as cm

pplt.rc.update({'xtick.major.size':6,
'xtick.minor.size':3.5,
'ytick.major.size':6,
'ytick.minor.size':3.5,
'tick.dir':'in',
'tick.labelsize':12,
'legend.columnspacing':0.1,
'legend.handletextpad':0.1,
'axes.labelsize':14,
'legend.fontsize':11,
#'cycle':'tab10'
})

图1

t=1e6*365*24*3600
c0=0.2e4
c1=4e4
h=200

d0=10
d1=-16
rsmow=0.0020052
c0D=(d0/1000+1)*c0*rsmow
c1D=(d1/1000+1)*c1*rsmow

fig,axs=pplt.subplots(ncols=2,nrows=1,figsize=(8,3),share=False)

#D0=np.array([10**2.24,10**1.11,np.exp(-7.3)/10000])[0]
#E=(np.array([374,334,118])*1e3)[0]
R=8.314
T=np.arange(200,601,1)
x=np.linspace(0,h,100)

x1,T1,j=np.meshgrid(x*1e-6,T+273.15,np.arange(0,300))
#D=D0*np.exp(-E/(R*T1))
DD=9e-9*np.exp(-49*4.184e3/R/T1) # diffusion coefficient of 2H
D=DD*np.sqrt(18/16) # diffusion coefficient of 1H

#######(a) dD varying with distance and temperature

sigma=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*D*t)
sumj=sigma.sum(axis=2)
c=c1+(c0-c1)*sumj

sigmaD=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*DD*t)
sumjD=sigmaD.sum(axis=2)
cD=c1D+(c0D-c1D)*sumjD

dD=(cD/c/rsmow-1)*1000

m=axs[0].contourf(x,T,dD,levels=np.linspace(-50,10,100),cmap='Oranges_r',extend='both')
axs[0].colorbar(m,loc='r',ticks=np.arange(-50,11,10),minorlocator=10,label=r'$\mathrm{\delta ^{18}O\ (‰)}$',width=0.15,pad=0.06,
                ticklabelsize=8,labelsize=10)
axs[0].text(0.97,0.9,'t=1 Ma',transform=axs[0].transAxes,fontsize=10,ha='right',c='k')
axs[0].format(xlim=(0,h),xlabel=r'$\mathrm{Distance\ ({\mu}m)}$',ylabel=r'$\mathrm{Temperature\ ({^o}C)}$')

##### (b)

t=10e6*365*24*3600

sigma=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*D*t)
sumj=sigma.sum(axis=2)
c=c1+(c0-c1)*sumj

sigmaD=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*DD*t)
sumjD=sigmaD.sum(axis=2)
cD=c1D+(c0D-c1D)*sumjD

dD=(cD/c/rsmow-1)*1000

m=axs[1].contourf(x,T,dD,levels=np.linspace(-50,10,100),cmap='Oranges_r',extend='both')
axs[1].colorbar(m,loc='r',ticks=np.arange(-50,11,10),minorlocator=10,label=r'$\mathrm{\delta ^{18}O\ (‰)}$',width=0.15,pad=0.06,
                ticklabelsize=8,labelsize=10)
axs[1].text(0.97,0.9,'t=10 Ma',transform=axs[1].transAxes,fontsize=10,ha='right',c='k')
axs[1].format(xlim=(0,h),xlabel=r'$\mathrm{Distance\ ({\mu}m)}$',ylabel=r'$\mathrm{Temperature\ ({^o}C)}$')

axs.format(abc='a',abcloc='ul',abc_kw={'fontsize':14})
#fig.savefig('d18O-diffusion.jpg',dpi=600)
   

图2

t=1e6*365*24*3600
c0=0.2e4
c1=4e4
h=200

d0=-50
d1=-60
rsmow=155.76e-6
c0D=(d0/1000+1)*c0*rsmow
c1D=(d1/1000+1)*c1*rsmow

fig,axs=pplt.subplots(ncols=2,nrows=1,figsize=(8,3),share=False)

#D0=np.array([10**2.24,10**1.11,np.exp(-7.3)/10000])[0]
#E=(np.array([374,334,118])*1e3)[0]
R=8.314
T=np.arange(0,201,1)
x=np.linspace(0,h,100)

x1,T1,j=np.meshgrid(x*1e-6,T+273.15,np.arange(0,300))
#D=D0*np.exp(-E/(R*T1))
DD=6.71e-13*np.exp(-80.5e3/R/T1) # diffusion coefficient of 2H
D=DD*np.sqrt(2) # diffusion coefficient of 1H

#######(a) dD varying with distance and temperature

sigma=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*D*t)
sumj=sigma.sum(axis=2)
c=c1+(c0-c1)*sumj

sigmaD=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*DD*t)
sumjD=sigmaD.sum(axis=2)
cD=c1D+(c0D-c1D)*sumjD

dD=(cD/c/rsmow-1)*1000

m=axs[0].contourf(x,T,dD,levels=np.linspace(-300,-50,100),cmap='Oranges_r',extend='both')
axs[0].colorbar(m,loc='r',ticks=np.arange(-300,-49,50),minorlocator=10,label=r'$\mathrm{\delta D\ (‰)}$',width=0.15,pad=0.06,
                ticklabelsize=8,labelsize=10)
axs[0].text(0.97,0.9,'t=1 Ma',transform=axs[0].transAxes,fontsize=10,ha='right',c='k')
axs[0].format(xlim=(0,h),xlabel=r'$\mathrm{Distance\ ({\mu}m)}$',ylabel=r'$\mathrm{Temperature\ ({^o}C)}$')

##### (b)

t=10e6*365*24*3600

sigma=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*D*t)
sumj=sigma.sum(axis=2)
c=c1+(c0-c1)*sumj

sigmaD=4/(2*j+1)/np.pi*np.sin(((2*j+1)*np.pi*x1)/h*1e6)*np.exp(-(((2*j+1)*np.pi)/h*1e6)**2*DD*t)
sumjD=sigmaD.sum(axis=2)
cD=c1D+(c0D-c1D)*sumjD

dD=(cD/c/rsmow-1)*1000

m=axs[1].contourf(x,T,dD,levels=np.linspace(-300,-50,100),cmap='Oranges_r',extend='both')
axs[1].colorbar(m,loc='r',ticks=np.arange(-300,-49,50),minorlocator=10,label=r'$\mathrm{\delta D\ (‰)}$',width=0.15,pad=0.06,
                ticklabelsize=8,labelsize=10)
axs[1].text(0.97,0.9,'t=10 Ma',transform=axs[1].transAxes,fontsize=10,ha='right',c='k')
axs[1].format(xlim=(0,h),xlabel=r'$\mathrm{Distance\ ({\mu}m)}$',ylabel=r'$\mathrm{Temperature\ ({^o}C)}$')

axs.format(abc='a',abcloc='ul',abc_kw={'fontsize':14})
#fig.savefig('dD-diffusion.jpg',dpi=600)