Minimal ray sphere intersection in C

sphere.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <png.h>
#include <assert.h>

struct vec3D {
   double x,y,z;
};

int width, height;
static double dot_product(const struct vec3D a, const struct vec3D b){
   return a.x*b.x+a.y*b.y+a.z*b.z;
}
static double magnitude(const struct vec3D a){
   return sqrt(dot_product(a,a));
}

int main(int argc, char *argv[]){
   double Pz, c, r;
   if (argc != 4) {
      printf("usage: %s Pz c r\n\n"
              "Pz = image plane z location\n"
              "c = sphere z center\n"
              "r = sphere radius\nusing default: Pz = 40, c = 202, r = 200\n" , argv[0]);
      Pz = 40.0; c = 202.0; r = 200.0;
   } else {Pz = atof(argv[1]), c = atof(argv[2]), r = atof(argv[3]);}
   width = 1000, height = 1000;

   FILE *fp = fopen("t.png", "wb");
   assert(fp!=NULL);

   png_structp png = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
   assert(png);
   png_infop info = png_create_info_struct(png);
   assert(info);
   int tem = setjmp(png_jmpbuf(png));
   assert(!tem);

   png_init_io(png, fp);
   png_set_IHDR(
      png,
      info,
      width, height,
      8,
      PNG_COLOR_TYPE_RGBA,
      PNG_INTERLACE_NONE,
      PNG_COMPRESSION_TYPE_DEFAULT,
      PNG_FILTER_TYPE_DEFAULT);
   png_write_info(png, info);

   png_bytep *row = NULL;
   row = malloc(sizeof(png_bytep)*height);
   assert(row!=NULL);

   for (int y = 0; y < height; y++){
      row[y] = malloc(sizeof(png_bytep)*width);
      assert(row[y]!=NULL);
   }

   struct vec3D light_source = {.x=-0.95,.y=0.0,.z=-0.31224}; /* magnitude = 1 */
   for (int y = 0; y < height; y++){
      png_bytep line = row[y];
      for (int x = 0; x < width; x++){
         double Px = x-(width/2), Py = (height/2)-y;
         struct vec3D ray = {Px,Py,Pz};
         double _magnitude = magnitude(ray);
         struct vec3D normalized_ray = {ray.x/_magnitude, ray.y/_magnitude, ray.z/_magnitude};

         struct vec3D sphere_center = {0.0,0.0,c};
         double sphere_radius = r;
         double v = dot_product(sphere_center, normalized_ray);
         double delta = sphere_radius*sphere_radius-((dot_product(sphere_center,sphere_center)-v*v));
         if (delta>0.0||delta==0.0){
            double d = sqrt(delta);
            struct vec3D intersection = {(v-d)*normalized_ray.x, (v-d)*normalized_ray.y, (v-d)*normalized_ray.z};
            /* center the temp_normal at origin */
            struct vec3D temp_normal = intersection;
            temp_normal.x -= sphere_center.x; temp_normal.y -= sphere_center.y;
            temp_normal.z -= sphere_center.z;

            /* normalized temp_normal */
            struct vec3D normal = {temp_normal.x/magnitude(temp_normal),
                                   temp_normal.y/magnitude(temp_normal),
                                   temp_normal.z/magnitude(temp_normal)};

            double _dot_product = dot_product(normal, light_source);
            if (_dot_product<0) _dot_product = 0.0;
            *(line+(x*4)+0)=(unsigned char)(_dot_product*255.0);
            *(line+(x*4)+1)=(unsigned char)(_dot_product*255.0);
            *(line+(x*4)+2)=(unsigned char)(_dot_product*255.0);
            *(line+(x*4)+3)=0xFF;

         }
         else if (delta<0.0){/* no intersection. AKA an imagninary sphere */
            *(line+(x*4)+0)=0x00;
            *(line+(x*4)+1)=0x00;
            *(line+(x*4)+2)=0x00;
            *(line+(x*4)+3)=0xFF;}
      }
   }

   png_write_image(png, row);
   png_write_end(png, NULL);

   for (int y = 0; y < height; y++) free(row[y]);
   free(row);

   fclose(fp);
   png_destroy_write_struct(&png, &info);
   return 0;
}