please covert this program from c++ to matlab

#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; } }

[KSSV]

Do you think people here will convert c++ code to MATLAB? Why you think people will spend time on your code? Also code looks lengthy; one needs to invest his time.

[Shailendra Singh Shah]

I need it very urgently. I know it's lengthy if someone is ready i will pay him and i will also stand with him while converting this. Thank you

[Rik]

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.