#include <math.h>
#include "/usr/local/pvm3/include/pvm3.h"

const double PI = 3.14159265358979323846;
const double TOL = 1E-6;

double timestep = 0.001;

// Projectile and Target Settings
double xtarget=800, rtarget = 20;
double rball = 0.1;
double aball = PI*rball*rball;
double mball = 1E2;

// Important Constants
const double g = 9.81, CD = 0.5, RHO=1.225;
double cd_rho_aball_over_2 = CD*RHO*aball/2.0;

// Math stuff and a place to store the answers before I send them back.
double project(double v0, double theta0);
void normalize(double *v);
double absval(double a);
double data[1000][2];

int main()
{
  int mytid, nproc, master, msgtype;

  int ppc, pcount=0;
  double vel, dtheta, hit, color;
  double theta = 0;

  mytid = pvm_mytid();

  // Find out which v0 we're working on and the number of points to compute.
  msgtype = 1729;
  pvm_recv(-1, msgtype);
  pvm_upkdouble(&vel, 1, 1);
  pvm_upkdouble(&dtheta, 1, 1);
  pvm_upkint(&ppc, 1, 1);

  // Go through the entire range, find the result of porject()-ing it,
  // and set a color value based on whether or not it hits the target.
  for (theta = 0; pcount < ppc; theta += dtheta, pcount++)
  {
    hit = project(vel, theta);

    data[pcount][0] = theta;
    if (hit >= 0.0)
    {
      color = ((1-(hit/rtarget))*0.75)+0.25;
      data[pcount][1] = color;
    }
    else
    {
      color = (1000+hit)/1000;
      data[pcount][1] = -1*color;
    }
  }

  // Send all the results back to the master program.
  msgtype = 44944;
  pvm_initsend(PvmDataDefault);
  pvm_pkdouble(&vel, 1, 1);
  pvm_pkdouble(&data[0][0], 2000, 1);
  master = pvm_parent();  
  pvm_send(master, msgtype);

  return 0;
}


double project(double v0, double theta0)
{
  double xp0=0.0, yp0=0.0;
  double xpos=0.0, ypos=0.0;
  double vlast[2], vnew[2];
  double force[2], d[2], dmag2, fmag, vmag, vang, i;
  double DT=timestep;

    // Calculate the trajectory.
    xpos=xp0, ypos=yp0, vnew[0]=v0*cos(theta0), vnew[1]=v0*sin(theta0);
    do {
      vlast[0]=vnew[0], vlast[1]=vnew[1];

      // Update position based on velocity
      if (ypos+vlast[1]*DT < 0)
      {
	xpos += (vlast[0] * absval(ypos/vlast[1]));
	ypos = 0;
      }
      else
	xpos += vlast[0]*DT, ypos += vlast[1]*DT;
      
      // GRAVITY - EARTH
      vnew[1] -= g*DT;

      // DRAG FORCE
      vmag = hypot(vlast[0], vlast[1]);
      vang = atan2(vlast[1], vlast[0]);
      fmag = cd_rho_aball_over_2*(vmag*vmag);
      vnew[0] -= (fmag/mball)*DT*cos(vang);
      vnew[1] -= (fmag/mball)*DT*sin(vang);

    } while (ypos > 0);
    if (absval(xpos-xtarget) <= rtarget)
      return absval(xpos-xtarget);
    return -1*(absval(xpos-xtarget));
}

void normalize(double *v)
{
  // Find the unit vector in the same direction as v.
  double dist = hypot(v[0], v[1]);
  v[0] /= dist;
  v[1] /= dist;
}

double absval(double a)
{
  // The abs() in stdlib.h doesn't do ints. ~shrug~
  if (a >= 0) return a;
  return (-1*a);
}

---

	My Supercomp front page.
	This page was created by Gary Sivek for 7th Period Supercomp.