#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cerrno>
#include <cctype>
#include <cmath>

static float to_radians(float d)
{
     return d / 180 * 3.141596;
}

char *getstr(FILE *fp)
{
     static char str[256];
     char *s = str;
     char c;

     while(1) {
	c = fgetc(fp);
	if(c == ';' || c == '\n')
	     break;
	*s++ = c;
     }
     *s = '\0';

     return str;
}

const float aspect = 1.0;

const float NO_t = -1.0;

const int NONE = 0;
const int SHAPE = 2;

typedef float real_t;

struct vec3f { /* Point in three dimensions... */
    real_t X, Y, Z;
    vec3f& operator=(const vec3f& p) {
        X = p.X;
        Y = p.Y;
        Z = p.Z;
        return *this;
    }
    vec3f& set(real_t X_, real_t Y_, real_t Z_)
    {
        X = X_;
        Y = Y_;
        Z = Z_;
        return *this;
    }
};

struct shape {
    vec3f C; /* Center of sphere */
    vec3f color;
    real_t R; /* (radius) */
};

struct vpoint { /* Viewpoint specification. */
    vec3f V; /* Origin of rays */
    real_t Yaw, Pit, Roll; /* Also known as Heading, Pitch, and Bank. */
    real_t VAng; /* Entire View angle, left to right. */
};

struct world { /* World specification */
    char *Title;
    shape *Shapes;
    vec3f BC; /* Background color */
    vpoint VP; /* Viewpoint */
    int NShapes; /* Number of polygons and shapes, */
    int XS, YS; /* Yaw and Pitch divisions per pixel */
    int df, lace;

    world() :
        Title(NULL),
        Shapes(NULL)
    {};
};

struct ray { /* Ray in 3D, parametric parameters. */
    vec3f O; /* Origin of the ray */
    vec3f d; /* dx, dy, and dz with respect to 't' */
};

vec3f vec3f_subtract(const vec3f& v1, const vec3f& v2)
{
    vec3f v;
    v.set(v1.X - v2.X, v1.Y - v2.Y, v1.Z - v2.Z);
    return v;
}

vec3f vec3f_add(const vec3f& v1, const vec3f& v2)
{
    vec3f v;
    v.set(v1.X + v2.X, v1.Y + v2.Y, v1.Z + v2.Z);
    return v;
}

vec3f vec3f_divide(const vec3f& v, real_t d)
{
    vec3f r;
    r.set(v.X / d, v.Y / d, v.Z / d);
    return r;
}

vec3f vec3f_multiply(const vec3f& v, real_t d)
{
    vec3f r;
    r.set(v.X * d, v.Y * d, v.Z * d);
    return r;
}

real_t vec3f_dot(const vec3f& v1, const vec3f& v2)
{
    return v1.X * v2.X + v1.Y * v2.Y + v1.Z * v2.Z;
}

/* Does ray intersect shape and where*/
bool shape_intersect(const shape& shp, const ray& ray, real_t *t)
{
    vec3f diff = vec3f_subtract(ray.O, shp.C);

    real_t b = 2 * vec3f_dot(ray.d, diff);

    real_t sqroot = b * b - 4 * (vec3f_dot(diff, diff) - shp.R * shp.R);
    if(sqroot < 0)
        return false;

    real_t t2 = (-b + sqrtf(sqroot)) / 2;

    if(t2 < 0)
         return false;

    real_t t1 = (-b - sqrtf(sqroot)) / 2;

    if(t1 < 0)
        *t = t2;
    else
        *t = t1;
    return true;
}

void FreeWorld(world *w)
{
    if(w->Shapes)
        delete[] w->Shapes;
    if(w->Title)
        free(w->Title);
    free(w);
}

int getint(FILE *fp, int *num)                  /* Call myio routine getstr() and */
                                                /*  convert result to integer */
{
    char *str;                                  /* String receiving from file */
 
    if((str = getstr(fp)) == NULL)
        return 0;
    *num = atoi(str);
    return 1;
}

real_t fltget(FILE *fp, real_t *num)            /* Call getstr() and convert result to*/
                                                /*  double-floating point */
{
    char *str;                                  /* Receipt from file */
 
    if((str = getstr(fp)) == NULL)
        return 0;
    *num = atof(str);
    return 1;
}

void usage() /* User-friendly help routine */
{
    fprintf(stderr, "usage: sphray.new <inputfilename>\n");
    fprintf(stderr, "A PPM file is written as output.");
}

vec3f Pixel(const ray& Ray, const world& Wd)
{
    int theone, thenum;
    real_t closest;
    vec3f color;

    theone = NONE;
    closest = 10000000.0;
    for(int cntr = 0; cntr < Wd.NShapes; cntr ++) {
        real_t distance;
        if(shape_intersect(Wd.Shapes[cntr], Ray, &distance)) {
            if((theone == NONE) || ((theone != NONE) && (distance < closest))) {
                theone = SHAPE;
                thenum = cntr;
                closest = distance;
            }
        }
    }
    if(theone == NONE)
        return Wd.BC;

    shape& tshp = Wd.Shapes[thenum];
    float brightness;
    if(Wd.df > 0)
        brightness = (tshp.C.Z - (Ray.O.Z + closest * Ray.d.Z)) / tshp.R + (random() % Wd.df - Wd.df / 2) / 100.0;
    else
        brightness = (tshp.C.Z - (Ray.O.Z + closest * Ray.d.Z)) / tshp.R;
    if(brightness < 0)
        brightness = 0;
    if(brightness > 1)
        brightness = 1;
    color = vec3f_multiply(tshp.color, brightness);
    return color;
}

void TraceImage(int width, int height, vec3f *image, const world& Wd)
{
    real_t indgx, indgy, RCOS, RSIN;

    RCOS = cos(Wd.VP.Roll) * cos(Wd.VP.Roll);
    RSIN = sin(Wd.VP.Roll) * sin(Wd.VP.Roll);

    indgx = Wd.VP.VAng / (width * Wd.XS);
    indgy = Wd.VP.VAng / (height * Wd.YS);

#pragma omp parallel for
    for(int yloop = 0; yloop < height; yloop++) {
        for(int xloop = 0; xloop < width; xloop++) {
            int cols = 0;
            vec3f color;
            color.set(0, 0, 0);
            for(int qkx = 0; qkx < Wd.XS; qkx++) {
                ray Ray; /* Ray that is being traced */
                real_t RawYaw, RawPit, YawAng, PitAng, PitCOS;
                RawYaw = Wd.VP.VAng / 2 - (xloop * Wd.XS + qkx) * indgx;
                for(int qky = 0; qky < Wd.YS; qky++) {
                    Ray.O = Wd.VP.V;
                    RawPit = Wd.VP.VAng / 2 - (yloop * Wd.YS + qky) * indgy;
                    YawAng = Wd.VP.Yaw + RawYaw * RCOS + RawPit * RSIN;
                    PitAng = Wd.VP.Pit + aspect * (-RawYaw * RSIN+RawPit * RCOS);
                    PitCOS = cos(PitAng);
                    Ray.d.set(sin(YawAng) * PitCOS, sin(PitAng), cos(YawAng) * PitCOS);
                    vec3f sample = Pixel(Ray, Wd);
                    color = vec3f_add(color, sample);
                    ++cols;
                }
            }
            image[yloop * width + xloop] = vec3f_divide(color, cols);
        }
        fprintf(stderr, "."); fflush(stderr);
    }
    fprintf(stderr, "\n");
}

world *GetWorld(char *fname)      /* Get world and return pointer. */
{
    FILE *fp;                     /* Pointer to file control structure. */
    char *inpstr;                 /* What I'm reading from the file */
    bool success;
                                                    /*  routine. */
    world *Wd;                         /* The World. */
   
    if((Wd = new world) == NULL) {
        fprintf(stderr, "Out of Mem; cannot allocate world structure.\n");
        return NULL;   
    }
    if((fp = fopen(fname, "r")) == NULL) {
        fprintf(stderr, "Cannot open file %s for input.\nE#%d\n", fname, errno);
        free(Wd); return NULL;
    }
    if((inpstr = getstr(fp)) == NULL) {
        fprintf(stderr, "Cannot read title.\n");
        fclose(fp);
        free(Wd);
        return NULL;
    }
    if((Wd->Title = (char *)malloc(strlen(inpstr) + 1)) == NULL) {
        fprintf(stderr, "No memory for world title.\n");
        goto nl2;
    }
    if(!strcmp(inpstr, ".")) {
        // FreeMem(Wd, sizeof(*Wd)); // XXX Why did I free here, and then set Wd->Title to NULL on the next line?  Plus continue reading!?
        Wd->Title = NULL;
    }
    else 
        strcpy(Wd->Title, inpstr);
    if(!getint(fp, &Wd->lace)) {
        fprintf(stderr, "*!LACE\n");
        goto nl2;
    }

    int wdth, lnth; // now ignored
    if(!getint(fp, &wdth)) {
        fprintf(stderr, "*!WDTH\n");
        goto nl2;
    }
    if(!getint(fp, &lnth)) {
        fprintf(stderr, "*!LNTH\n");
        goto nl2;
    }

    if(!getint(fp, &Wd->XS)) {
        fprintf(stderr, "*!XS\n");
        goto nl2;
    }
    if(!getint(fp, &Wd->YS)) {
        fprintf(stderr, "*!YS\n");
        goto nl2;
    }
    int r, g, b, brt;
    if(!getint(fp, &r)) {
        fprintf(stderr, "*!BCr\n");
        goto nl2;
    }
    if(!getint(fp, &g)) {
        fprintf(stderr, "*!BCg\n");
        goto nl2;
    }
    if(!getint(fp, &b)) {
        fprintf(stderr, "*!BCb\n");
        goto nl2;
    }
    if(!getint(fp, &brt)) {
        fprintf(stderr, "*!BCb\n");
        goto nl2;
    }
    Wd->BC.X = r / 16.0 * brt / 100.0;
    Wd->BC.Y = g / 16.0 * brt / 100.0;
    Wd->BC.Z = b / 16.0 * brt / 100.0;
    if(!getint(fp, &Wd->df)) {
        fprintf(stderr, "*!DIFFUSION\n");
        goto nl2;
    }
    if(!getint(fp, &Wd->NShapes)) {
        fprintf(stderr, "*!#Sources\n");
        goto nl2;
    }
    if((Wd->Shapes = new shape[Wd->NShapes]) == NULL) {
        fprintf(stderr, "No memory for light shapes.\n");
        goto nl2;
    }
    success = true;
    for(int cntr = 0; cntr < Wd->NShapes; cntr++) {
        shape& tshape = Wd->Shapes[cntr];
        int r, g, b, brt;
        success = fltget(fp, &tshape.C.X) && fltget(fp, &tshape.C.Y);
        success = success && fltget(fp, &tshape.C.Z) && fltget(fp, &tshape.R);
        success = success && getint(fp, &r) && getint(fp, &g);
        success = success && getint(fp, &b) && getint(fp, &brt);
        if(!success) {
            fprintf(stderr, "*!Source #%d\n", cntr);
            goto nl4;
        }
        tshape.color.X = r / 16.0 * brt / 100.0;
        tshape.color.Y = g / 16.0 * brt / 100.0;
        tshape.color.Z = b / 16.0 * brt / 100.0;
    }
    success = success && fltget(fp, &Wd->VP.V.X) && fltget(fp, &Wd->VP.V.Y);
    success = success && fltget(fp, &Wd->VP.V.Z) && fltget(fp, &Wd->VP.Yaw);
    success = success && fltget(fp, &Wd->VP.Pit) && fltget(fp, &Wd->VP.Roll);
    success = success && fltget(fp, &Wd->VP.VAng);
    if(!success) {
        fprintf(stderr, "*!Viewpoint\n");
        goto nl4;
    }
    fclose(fp);

    Wd->VP.Pit = to_radians(Wd->VP.Pit);
    Wd->VP.Yaw = to_radians(Wd->VP.Yaw);
    Wd->VP.Roll = to_radians(Wd->VP.Roll);
    Wd->VP.VAng = to_radians(Wd->VP.VAng);

    return Wd;

nl4:
    free(Wd->Shapes);
nl2:
    if(Wd->Title != NULL)
        free(Wd->Title);
    fclose(fp);
    free(Wd);
    return NULL;
}

int main(int argc, char **argv)
{
    world *World;

    if(argc < 2) {
        usage();
        exit(EXIT_FAILURE);
    }

    if((!strcmp(argv[1], "-h")) || (!strcmp(argv[1], "--help"))) {
        usage();
        exit(EXIT_FAILURE);
    }

    if((World = (world *)GetWorld(argv[1])) == NULL) {
        fprintf(stderr, "Cannot set up world.\n");
        exit(EXIT_FAILURE);
    }

    int width = 512, height = 512;
    vec3f *image = new vec3f[width * height];

    TraceImage(width, height, image, *World);

    printf("P6 %d %d 255\n", width, height);
    for(int j = 0; j < height; j++)
         for(int i = 0; i < width; i++) {
              vec3f& c = image[j * width + i];
              printf("%c%c%c", (int)(c.X * 255), (int)(c.Y * 255), (int)(c.Z * 255));
         }

    delete[] image;
    FreeWorld(World);
    exit(EXIT_SUCCESS);
}