Néhány módosítás kell még, de ez már egész jó lett.

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


#define xh 700
#define yh 500
unsigned char bitmap[xh*yh*4];


void pont(int x,int y,int r,int g,int b)
{
    if(x<0) return;
    if(x>=xh) return;
    if(y<0) return;
    if(y>=yh) return;
    
    int address=(y*xh+x)*3;
    bitmap[address+2]=r;
    bitmap[address+1]=g;
    bitmap[address]=b;
}
struct vec3
{
    float x,y,z;
    
    vec3() {x=y=z=0.0;}
    vec3(float x3,float y3,float z3) {x=x3;y=y3;z=z3;}
    vec3 operator +(float s3) {vec3 v3;v3.x=x+s3;v3.y=y+s3;v3.z=z+s3;return v3;}
    vec3 operator -(float s3) {vec3 v3;v3.x=x-s3;v3.y=y-s3;v3.z=z-s3;return v3;}
    vec3 operator *(float s3) {vec3 v3;v3.x=x*s3;v3.y=y*s3;v3.z=z*s3;return v3;}
    vec3 operator /(float s3) {vec3 v3;v3.x=x/s3;v3.y=y/s3;v3.z=z/s3;return v3;}

    vec3 operator +(vec3 v4) {vec3 v3;v3.x=x+v4.x;v3.y=y+v4.y;v3.z=z+v4.z;return v3;}
    vec3 operator -(vec3 v4) {vec3 v3;v3.x=x-v4.x;v3.y=y-v4.y;v3.z=z-v4.z;return v3;}
    vec3 operator *(vec3 v4) {vec3 v3;v3.x=x*v4.x;v3.y=y*v4.y;v3.z=z*v4.z;return v3;}
    vec3 operator /(vec3 v4) {vec3 v3;v3.x=x/v4.x;v3.y=y/v4.y;v3.z=z/v4.z;return v3;}
};

float dot(vec3 v1,vec3 v2) {return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;}

vec3 cross(vec3 v1,vec3 v2)
{
    vec3 v3;
    v3.x=v1.y*v2.z - v1.z*v2.y;
    v3.y=v1.z*v2.x - v1.x*v2.z;
    v3.z=v1.x*v2.y - v1.y*v2.x;
    
    return v3;
}
vec3 normalize(vec3 v3)
{
    return v3/sqrt(dot(v3,v3));
}

struct mat3
{
    vec3 mx,my,mz,mw;
};



float radian=M_PI/180.0;
mat3 view;


void vonal3d(vec3 v1,vec3 v2,int r,int g,int b)
{
   for(int i=0;i<1000;i++)
    {
        vec3 v4,v3=v1+(v2-v1)*(float)i/1000.0;

        v4.x=dot(v3,view.mx);
        v4.y=dot(v3,view.my);
        v4.z=dot(v3,view.mz);
        v4=v4+view.mw;

        pont(xh/2+(int)v4.x,yh/2+(int)v4.y, r,g,b);
    }
}
void set_view_polarKR(float fi1,float fi2,float skala,vec3 offset)
{
    fi1*=radian;
    fi2*=radian;
    
    view.my=vec3(0,1,0);
    view.mz=vec3(cos(fi1)*cos(fi2), sin(fi2), sin(fi1)*cos(fi2));
    view.mx=cross(view.mz,view.my);
    view.my=cross(view.mx,view.mz);

    view.mx=normalize(view.mx);
    view.my=normalize(view.my);
    view.mz=normalize(view.mz);
    
    view.mx=view.mx*skala;
    view.my=view.my*skala;
    view.mz=view.mz*skala;
    view.mw=offset;
    
}
int main()
{
    set_view_polarKR(40,30,0.025,vec3(0,-250,0));
    
    

    float c=1;
//    float v=0.999973*c;//gamma~137
    float v=0.8*c;
    float gamma=1/sqrt(1-v*v/(c*c));
    float x1,t1,x2,t2;
    
    for(int ii=0;ii<137*8;ii++)    
    {
        float fi3=ii*M_PI*2/137.0;
        
        float ti=(float)ii*8.0;
        vec3 p0=vec3(0,ti,0);
        vec3 p1=vec3(500.0*cos(fi3),500+ti, 500.0*sin(fi3));
        vec3 p2,p02;

        
        x1=p1.x;
        t1=p1.y;//ido
        x2=(x1-v*t1)*gamma;
        t2=(t1-v*x1/(c*c))*gamma;
        p2.x=x2;
        p2.y=t2;//ido
        p2.z=p1.z;
        
        x1=p0.x;
        t1=p0.y;//ido
        x2=(x1-v*t1)*gamma;
        t2=(t1-v*x1/(c*c))*gamma;
        p02.x=x2;
        p02.y=t2;//ido
        p02.z=p0.z;
        
        vonal3d(p0,p1,  255-ii,0,ii);
        vonal3d(p02,p2,  255,255-ii,0);
    }


    
    short fejlec[]={0,2,0,0,0,0,xh,yh,3*8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    FILE *file=fopen("save.tga","wb");
    fwrite(fejlec,1,18,file);
    fwrite(bitmap,1,xh*yh*3,file);
    fclose(file);

    return 0;
}