MATLAB Answers

Info

この質問は閉じられています。 編集または回答するには再度開いてください。

please covert this program from c++ to matlab

3 ビュー (過去 30 日間)
Shailendra Singh Shah
Shailendra Singh Shah 2021 年 3 月 1 日
終了済み: Rik 2021 年 3 月 1 日
この 質問 は 2 人のコントリビューターによってフラグが設定されました
#include<stdio.h> #include<math.h> #include<iostream.h> #include<stdlib.h> #include<fstream.h>
double tetar_sand,tetas_sand,ks_sand,beta_sand,gamma_sand,alfa_sand,pa_sand; double tetar_clay,tetas_clay,ks_clay,beta_clay,gamma_clay,alfa_clay,pa_clay; double tetar_van,tetas_van,ks_van,beta_van,gamma_van,alfa_van,pa_van; double qwin,qwou,vwp,ratio,baler,hs;
void main() { double fun_kt(double,int); double fun_kh(double,int); double fun_ct(double,int); double fun_ht(double,int); double fun_th(double,int); int nnode, n1,mat_type,kblock,kboun,nlim,nsteps,itermx,iflag,iter; double cc,dz,dt,ft,fto,dz2,rate,hbn,hb1,bet,ddiv,dmul,dtmax,zmin,eps,tfinal,tiniz; double tprint,dtprint,dtmin,vwi,cwin,cwou,dto,time1,dhmax,uh; double h[500],fk[500],km[500],a[500],b[500],c[500],r[500],u[500],ho[500],gam[500],z[500 ]; // variables defination ks_sand=9.22E-3; ks_clay=1.23E-5; ks_van=8.67E-5; pa_sand=1.175E+6; pa_clay=124.6; pa_van=0.0; alfa_sand=0.0335; alfa_clay=739.0; alfa_van=8.50E-2; beta_sand=2.0; beta_clay=1.77; beta_van=1.246475; gamma_sand=0.5; gamma_clay=4.0; gamma_van=0.12; tetas_sand=0.381; tetas_clay=0.495; tetas_van=0.5315; tetar_sand=0.102; tetar_clay=0.124; tetar_van=0.05325; mat_type=3; nnode=31; kblock=1; kboun=1;
dtmin=1.e-5; dtprint=3600.0; tprint=3600.0; itermx=30; nsteps=12000; tiniz=0.0; tfinal=4000.0; eps=1.e-8; zmin=0.0; dz=2.0; dt=1.0e-6; dtmax=10; // cahnged it from 10 to dmul=1.1; ddiv=0.5; nlim=10; ofstream fout, jout; fout.open ("RichaM2CP_out_5densgrid.xls"); jout.open ("RichaM2CP1_5_out.xls"); cout<<"PLEASE WAIT PROGRAMM IS RUNNING"<<endl; jout.width(10); jout<<"steps no"; jout.width(20); jout<<"Time"; jout.width(15); jout<<"Iteration"; jout.width(20); jout<<"DT"; jout.width(20); jout<<"DHMAX"<<"\n\n"; z[1]=zmin; for (int i=2; i <= nnode; i++) { z[i]=z[i-1]+dz; }
ho[1]=-4600; for (i=2; i <= nnode; i++) { ho[i]=-4600; } vwi=0.0; cwin=0.0; cwou=0.0; for (i=1; i <= nnode; i++) { vwi=vwi+dz*fun_th(ho[i],mat_type); } if(kboun==1)rate=4.62e-5; fout<<"\n\nINITIAL CONDITION TABLE\n\n"<<endl; fout.width(5); fout<<"Node"; fout.width(8); fout<<"Depth"; fout.width(15); fout<<"Pressure Head"; fout.width(25); fout<<"Moisture Content"; fout.width(25); fout<<"Conductivity"<<"\n"<<endl;
for (i=1; i <= nnode; i++) { h[i]=ho[i]; fout.width(5); fout<<i; fout.width(8); fout<<z[i]; fout.width(15); fout<<ho[i]; fout.width(25); fout<<fun_th(ho[i],mat_type); fout.width(25); fout<<fun_kh(ho[i],mat_type)<<"\n\n\n"; } fout<<"Initial volume of water: "<<vwi<<endl; if(kboun==1)cout<<"Prefixed rate: "<<rate<<endl; hb1=0.0; hbn=0.0; n1=nnode-1; time1=tiniz; dto=dt; dz2=dz*dz; for (int n=1; n <= nsteps; n++) { // open for loop iter=1; iflag=0; ASSFD:
for (i=1; i <= nnode; i++) { fk[i]=fun_kh(h[i],mat_type); } for (i=1; i <= n1; i++) { if(kblock==1) km[i]=0.5*(fk[i]+fk[i+1]); else if(kblock==2) km[i]=2.0*(fk[i]*fk[i+1])/(fk[i]+fk[i+1]); else if(kblock==3) km[i]=sqrt(fk[i]*fk[i+1]); else if(kblock==4) { if(h[i]>=h[i+1]) km[i]=fk[i]; else if(h[i]<h[i+1]) km[i]=fk[i+1]; } } for (i=2; i <= n1; i++) { cc=fun_ct(h[i],mat_type); a[i]=-km[i-1]/dz2; b[i]=cc/dt+(km[i-1]+km[i])/dz2; c[i]=-km[i]/dz2; r[i]=(km[i-1]*(h[i-1]-h[i])+km[i]*(h[i+1]-h[i]))/dz2-(km[i]km[i-1])/dz; ft=fun_th(h[i],mat_type); fto=fun_th(ho[i],mat_type); r[i]=r[i]-(ft-fto)/dt; }
// boundary conditions if(kboun==0) { r[1]=hb1; r[2]=r[2]-a[2]*hb1; r[n1]=r[n1]-c[n1]*hbn; r[nnode]=hbn; a[2]=0.0; a[nnode]=0.0; b[1]=1.0; b[nnode]=1.0; c[1]=0.0; c[n1]=0.0; } else if(kboun==1) { r[nnode]=hbn; b[nnode]=1.0; a[nnode]=0.0; r[n1]=r[n1]-c[n1]*hbn; c[n1]=0.0; cc=fun_ct(h[1],mat_type); b[1]=cc/dt+km[1]/dz2; c[1]=-km[1]/dz2; ft=fun_th(h[1],mat_type); fto=fun_th(ho[1],mat_type); r[1]=km[1]*(h[2]-h[1])/dz2-km[1]/dz-(ft-fto)/dt+rate/dz; }
// tridiagonal matrix solution if(b[1]==0.0) { cout<<"divide by zero error"<<endl; exit(0); } bet=b[1]; u[1]=r[1]/bet; for (int j=2; j <= nnode; j++) { gam[j]=c[j-1]/bet; bet=b[j]-a[j]*gam[j]; if(bet==0.0)exit(0); u[j]=(r[j]-a[j]*u[j-1])/bet; } for (j=nnode-1; j >= 1; j--) { u[j]=u[j]-gam[j+1]*u[j+1]; } //close tridiagonal matrix solution
dhmax=1.e-30; for (i=1; i <= nnode; i++) { uh=u[i]/h[i]; if(fabs(uh)>dhmax)dhmax=fabs(uh); } if(dhmax > eps && iter<itermx) { for (i=1; i <= nnode; i++) { h[i]=h[i]+u[i];
} iter=iter+1; goto ASSFD; } else if(dhmax>eps && iter>=itermx && dt>=dtmin) { dt=dt*ddiv; iflag=iflag+1; for (i=1; i <= nnode; i++) { h[i]=h[i]+u[i]; h[i]=ho[i]+(h[i]-ho[i])*dt/dto; h[i]=ho[i]; } iter=1; goto ASSFD; } else if(dhmax>eps && iter>=itermx && dt<=dtmin) { cout<<"convergence not reached"<<endl; exit(0); } else if(dhmax<=eps) { time1=time1+dt; jout.width(10); jout<<n; jout.width(20); jout<<time1; jout.width(15); jout<<iter; jout.width(20); jout<<dt; jout.width(20); jout<<dhmax<<"\n\n"; if(time1>=tprint) { fout<<"\nReport time: (hr) "<<time1/3600<<endl; for (i=1; i <= nnode; i++) { h[i]=h[i]+u[i]; } } n1=nnode-1; qwin=km[1]*((h[1]-h[2])/dz+1.0); qwou=km[n1]*((h[n1]-h[nnode])/dz+1.0); cwin=cwin+qwin*dt; cwou=cwou+qwou*dt; if(time1>=tprint) { vwp=0.0; for (i=1; i <= nnode; i++) { vwp=vwp+fun_th(h[i],mat_type)*dz; } ratio=(vwp-vwi)/(cwin-cwou); baler=100.0*(1.0-ratio); fout<<"\nWater entered in step: "<<qwin<<endl; fout<<"\nCumulative water entered: "<<cwin<<endl; fout<<"\nWater out in the step: "<<qwou<<endl;
fout<<"\nCumulative water out: "<<cwou<<endl; fout<<"\nInitial water volume: "<<vwi<<endl; fout<<"\nActual water volume: "<<vwp<<endl; fout<<"\nWater balance error(%): "<<baler<<endl; fout<<"\nNumber of time steps: "<<n<<"\n\n"<<endl; fout.width(5); fout<<"Node"; fout.width(12); fout<<"Soil Type"; fout.width(12); fout<<"Depth"; fout.width(25); fout<<"Pressure Head"; fout.width(25); fout<<"Moisture Content"<<"\n"; for (i=1; i <= nnode; i++) { fout.width(5); fout<<i; fout.width(12); fout<<mat_type; fout.width(12); fout<<z[i]; fout.width(25); fout<<h[i]; fout.width(25); fout<<fun_th(h[i],mat_type)<<"\n"; } tprint=tprint+dtprint; } dto=dt; if(iter<=nlim && iflag==0) dt=dt*dmul; if(dt>dtmax) dt=dtmax; if((time1+dt) >tprint)dt=time1+dt-tprint; for (i=1; i <= nnode; i++) { hs=h[i]+(h[i]-ho[i])*dt/dto; ho[i]=h[i]; h[i]=hs; } }//close else if loop }// close for loop fout.close(); jout.close(); }// end of main
double fun_kt(double x,int mat_type) { double kt,s_sand,r_sand,a_sand,d_sand; double s_clay,r_clay,e_clay,d_clay; double n_van,m_van,um_van,s_van,a_van,aq_van; switch(mat_type) { case 1: if(x<=tetar_sand) { kt=0.0; return(kt); } else if(x >= tetas_sand) { kt=ks_sand; return(kt); } s_sand = (x-tetar_sand)/(tetas_sand-tetar_sand); r_sand = beta_sand/gamma_sand; a_sand = alfa_sand/(s_sand-alfa_sand); d_sand = pa_sand+pow(a_sand,r_sand); kt=ks_sand*pa_sand/d_sand; return(kt); break; case 2: if(x<=tetar_clay) { kt=0.0; return(kt); } else if(x >= tetas_clay) { kt=ks_clay; return(kt); } s_clay = (x-tetar_clay)/(tetas_clay-tetar_clay); r_clay = beta_clay/gamma_clay; e_clay = pow((alfa_clay/s_clay),r_clay); d_clay = pa_clay + exp(e_clay); kt = ks_clay * pa_clay/d_clay; return(kt); break; case 3: if(x<=tetar_van) { kt=0.0; return(kt); } else if(x >= tetas_van) { kt=ks_van; return(kt); } n_van=beta_van; m_van=1.0-1.0/n_van; um_van=1.0/m_van; s_van=(x-tetar_van)/(tetas_van-tetar_van); a_van=1.0-pow((1.0-pow(s_van,um_van)),m_van);
aq_van=a_van*a_van; kt=ks_van*sqrt(s_van)*aq_van; return(kt); break; } }
double fun_ct(double x,int mat_type) { double ct,a_sand,gamma_sand1,b_sand; double a_clay,gamma_clay1,b1_clay,b2_clay,b_clay; double n_van,m_van,a1_van,a2_van,a_van,rn_van,rd_van; if((x>=0.0)||(x<-1.0e35)) { ct=0.0; return(ct); } switch(mat_type) { case 1: a_sand=pow((alfa_sand+pow(fabs(x),gamma_sand)),2); gamma_sand1=gamma_sand-1.0; b_sand=alfa_sand*(tetas_sandtetar_sand)*gamma_sand*pow(fabs(x),gamma_sand1); ct=b_sand/a_sand; return(ct); break; case 2: if(x>=-1.0) { ct=0.0; return(ct); break; } a_clay=pow((alfa_clay+log(fabs(x))*gamma_clay),2); gamma_clay1=gamma_clay-1.0; b1_clay=alfa_clay*(tetas_clay-tetar_clay)*gamma_clay; b2_clay=pow(log(fabs(x)),gamma_clay1)/fabs(x); b_clay=b1_clay*b2_clay; ct=b_clay/a_clay; return(ct); break; case 3: n_van=beta_van; m_van=1.0-1.0/n_van;
a_van=alfa_van*fabs(x); a1_van=pow(a_van,n_van); a2_van=pow(a_van,(n_van-1.0)); rn_van=n_van*m_van*(tetas_clay-tetar_clay)*a2_van*alfa_van; rd_van=pow((1.0+a1_van),(m_van+1.0)); ct=rn_van/rd_van; return(ct); break; } }
double fun_ht(double x,int mat_type) { double ht,s_sand,a_sand,s_clay,ug_clay,e_clay; double n_van,m_van,s_van,a_van; switch(mat_type) { case 1: if(x<=tetar_sand) { ht=-1.0e30; return(ht); } else if(x>=tetas_sand) { ht=0.0; return(ht); } s_sand = (x-tetar_sand)/(tetas_sand-tetar_sand); a_sand = alfa_sand/s_sand-alfa_sand; ht=-pow(a_sand,(1.0/gamma_sand)); return(ht); case 2: if(x<=tetar_clay) { ht=-1.0e30; return(ht); } else if(x>=tetas_clay) { ht=0.0; return(ht); } s_clay = (x-tetar_clay)/(tetas_clay-tetar_clay); ug_clay=1.0/gamma_clay;
e_clay=pow((alfa_clay/s_clay-alfa_clay),ug_clay); ht=-exp(e_clay); return(ht); case 3: if(x<=tetar_van) { ht=-1.0e30; return(ht); } else if(x>=tetas_van) { ht=0.0; return(ht); } n_van=beta_van; m_van=1.0-1.0/n_van; s_van = (x-tetar_van)/(tetas_van-tetar_van); a_van =pow(s_van,(-1.0/m_van))-1.0; ht=-pow(a_van,(1.0/n_van))/alfa_van; return(ht); } }
double fun_kh(double x,int mat_type) { double kh,n_van,m_van,b_van,rd_van,rn_van,alfn_van,alfn_van1; switch(mat_type) { case 1: if(x>=0.0) { kh=ks_sand; return(kh); } kh=ks_sand*pa_sand/(pa_sand+pow(fabs(x),beta_sand)); return(kh); case 2: if(x>=0.0) { kh=ks_van; return(kh); } kh=ks_clay*pa_clay/(pa_clay+pow(fabs(x),beta_clay)); return(kh);
case 3: if(x>=0.0) { kh=ks_van; return(kh); } n_van=beta_van; m_van=1.0-1.0/n_van; b_van=alfa_van*fabs(x); alfn_van=pow(b_van,n_van); alfn_van1=pow(b_van,(n_van-1)); rn_van=pow((1.0-alfn_van1*pow((1.0+alfn_van),-m_van)),2); rd_van=pow((1.0+alfn_van),(m_van/2)); kh=ks_van*rn_van/rd_van; return(kh); } }
double fun_th(double x,int mat_type) { double th,a_sand,b_sand,a_clay,b_clay; double n_van,m_van,rd_van,alfn_van; switch(mat_type) { case 1: if(x>=0.0) { th=tetas_sand; return(th); } else if(x<-1.0e4) { th=tetar_sand; return(th); } a_sand=alfa_sand*(tetas_sand-tetar_sand); b_sand=alfa_sand+pow(fabs(x),gamma_sand); th=a_sand/b_sand+tetar_sand; return(th); break; case 2: if(x>=0.0) { th=tetas_clay; return(th); } else if(x<-1.0e4) { th=tetar_clay; return(th); }
else if(x>=-1.0) { th=tetas_clay; return(th); } a_clay=alfa_clay*(tetas_clay-tetar_clay); b_clay=alfa_clay+pow(log(fabs(x)),gamma_sand); th=a_clay/b_clay+tetar_clay; return(th); break; case 3: if(x>=0.0) { th=tetas_van; return(th); } else if(x<-1.0e4) { th=tetar_van; return(th); } n_van=beta_van; m_van=1.0-1.0/n_van; alfn_van=pow((alfa_van*fabs(x)),n_van); rd_van=pow((1.0+alfn_van),m_van); th=(tetas_van-tetar_van)/rd_van+tetar_van; return(th); break; } }
  3 件のコメント
Rik
Rik 2021 年 3 月 1 日
You can hire a consultant if you like (Mathworks provides such services, and others do as well, they're a quick Google search away). If you want help on this forum, have a read here and here. It will greatly improve your chances of getting an answer.

回答 (0 件)

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by