MATLAB Answers

Info

This question is closed. Reopen it to edit or answer.

please covert this program from c++ to matlab

22 views (last 30 days)
Shailendra Singh Shah
Shailendra Singh Shah on 1 Mar 2021
Closed: Rik on 1 Mar 2021
This question was flagged by 2 contributors
#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 Comments
Rik
Rik on 1 Mar 2021
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.

Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!