Someone

float radon_func(float rho, float theta, struct einfo w){  float xp1, xm1, yp1, ct, st, tt, xm, xp, ym, yp;  float base, gamma, a, b, I1, I2, I3, I4, y1, y0;  rho -= w.x0 * cos(theta) + w.y0 * sin(theta);  theta -= w.phi;  a = w.a;  b = w.b;  while( theta >= PI ) {    rho = -rho;    theta -= PI;  }  while( theta < 0 ) {    rho = -rho;    theta += PI;  }  ct = cos(theta);  st = sin(theta);  gamma = hypot(a * ct, b * st);  rho = rho / gamma;  if( theta >= PI2 ) theta = atan(b / a * st / ct) + PI;  else theta = atan(b / a * st / ct);  switch( w.phtype ) {    case 1 : /* Circular disc */      if( (tt = 1.0 - rho * rho) > 0.0 ) base = 2.0 * sqrt(tt);      else base = 0.0;      break;    case 2 : /* Square */      if( fabs(rho) > sqrt2 ) base = 0.0;      else {        if( rho < 0 ) {          rho = -rho;          theta += PI;        }        while( theta >= 2 * PI )          theta -= 2 * PI;        if( theta >= PI ) theta = 2 * PI - theta;        if( theta >= PI * 0.5 ) theta = PI - theta;        if( theta >= PI * 0.25 ) theta = PI * 0.5 - theta; /* Now theta <PI/4 */        ct = cos(theta);        st = sin(theta);        tt = st / ct;        xp1 = rho / ct - tt;        xm1 = rho / ct + tt;        if( theta < 1e-20 ) if( rho < 1 ) base = 2;        else base = 0.0;        else {          yp1 = (rho - ct) / st;          if( xp1 > 1 ) base = 0.0;          else if( xm1 < 1 ) base = 2 / ct;          else base = hypot(1 - xp1, 1 - yp1);        }      }      break;    case 3 : /* Gaussian bell */      base = exp(-rho * rho) * sqrt(PI);      break;    case 4 : /* Even side triangle */      if( rho < 0 ) {        rho = -rho;        theta += M_PI;      }      if( rho > 1 ) base = 0.0;      else {        while( theta > topi3 )          theta -= topi3;        while( theta < 0.0 )          theta += topi3;        if( theta > pi3 ) theta = topi3 - theta;        ct = cos(theta);        st = sin(theta);        if( rho > ct ) base = 0.0;        else {          xp = (rho + st / sqrt3) / (st / sqrt3 + ct);          yp = (xp - 1.0) / sqrt3;          if( (rho + 0.5 * ct) > (st * sqrt3_2) ) {            xm = (rho - st / sqrt3) / (ct - st / sqrt3);            ym = -(xm - 1.0) / sqrt3;            base = hypot((xp - xm), (yp - ym));          } else base = hypot((xp + 0.5), (yp - (rho + 0.5 * ct) / st));        }      }      break;    case 5 : /* Noise is added */      base = gasdev(idum);      break;    case 6 : /* "Pyramid" */      if( fabs(rho) > sqrt2 ) base = 0.0;      else {        if( rho < 0 ) rho = -rho;        while( theta >= 2 * PI )          theta -= 2 * PI;        if( theta >= PI ) theta = 2 * PI - theta;        if( theta >= PI * 0.5 ) theta = PI - theta;        if( theta >= PI * 0.25 ) theta = PI * 0.5 - theta; /* Now theta <PI/4 */        ct = cos(theta);        st = sin(theta);        if( rho > st + ct ) base = 0.0;        else if( st < 1e-7 ) base = (1 - rho);        else {          if( rho > -st + ct ) /* y1 > -1 */          {            y1 = (rho - ct) / st;            I1 = ct * (1 - y1);            if( y1 < 0 ) I2 = -0.5 * ct * (1 + y1 * y1);            else I2 = -0.5 * ct * (1 - y1 * y1);            if( rho > st ) {              I3 = rho * (0.5 * y1 - 1) + 0.5 * (st + y1 * ct);              if( y1 < 0 ) I4 = 0.5 * rho + (y1 * y1 * (ct + 0.5 * rho) - st) / 3.0;              else I4 = 0.5 * rho - (y1 * y1 * (ct + 0.5 * rho) + st) / 3.0;            } else {              y0 = rho / st;              I3 = rho * (1 - y0 + 0.5 * y1) + 0.5 * (y1 * ct - st);              if( y1 < 0 ) I4 = (st + y0 * y0 * rho + y1 * y1 * (ct + 0.5 * rho)) / 3 - 0.5 * rho;              else I4 = (st + y0 * y0 * rho - y1 * y1 * (ct + 0.5 * rho)) / 3 - 0.5 * rho;            }          } else /* y1 < -1 */          {            I1 = 2 * ct;            I2 = -ct;            if( rho > st ) {              I3 = -2 * rho;              I4 = rho;            } else {              y0 = rho / st;              I3 = -st - rho * rho / st;              I4 = (2 * st + rho * y0 * y0) / 3;            }          }          base = (I1 + I2 + I3 + I4) / (ct * ct);        }      }      break;    default :      Print(_DNormal, "Function of wrong type detected\n");      exit(1);  }  return base * w.pow * a * b / gamma;}