Weird stuff in InitGasVelocities
At the end of that function, we resset vtheta
for (int i = 0; i < nr; i++) {
for (int h=0; h< ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
real rg = Rmed[i]*sin(phiMed[h]);
vt[l] -= *disk.referential().omega*rg;
}
}
}
So what's the point in editing it before that ?
if (viscosity.type == VT_ALPHA || viscosity.dead){
if(i>2 ){
real gradphi=(1./Rmed[i]-1./Rmed[i-1])/(Rmed[i]-Rmed[i-1]);
real gradp=(pres[l]-pres[l-ni*ns])/(Rmed[i]-Rmed[i-1])/(rho[l]+rho[l-ni*ns])*2.0; // this is not good for a log-grid
if((-gradphi+gradp) >=0.) {
real vtinf= sqrt((-gradphi+gradp)*ri);// *disk.referential().omega added later
vt[l]=2.0*vtinf-vt[l-ni*ns];
} else {
printf("In PframeForce warning problems in computing azimuthal velocity: %d\t%d\t%d\t%f\n",i,j,h,-gradphi+gradp);
}
} else if(i==2) {
real gradphi=(1./Rmed[i]-1./Rmed[i-1])/(Rmed[i]-Rmed[i-1]);
real gradp=(pres[l]-pres[l-ni*ns])/(Rmed[i]-Rmed[i-1])/(rho[l]+rho[l-ni*ns])*2.0; // this is not good for a log-grid
real vtinf= sqrt((-gradphi+gradp)*ri); //*disk.referential().omega added later
vt[l-ni*ns] = 2*vtinf/(sqrt(Rmed[1]/Rmed[2])+1);
vt[l] = vt[l-ni*ns]*sqrt(Rmed[1]/Rmed[2]);
vt[l-2*ni*ns] = vt[l-ni*ns]*sqrt(Rmed[1]/Rmed[0]);
}
} else {
vt[l] = omega*rg*
sqrt(1.0-(pow(aspectRatio[i],2.0)*
(2.+disk.physic().density.slope-disk.physic().flaringIndex)));
}