/* Steven Andrews, 8/25/02 */
/* See documentation called Plot3D doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "DiskIO.h"
#include "Plot.h"
#include "Plot3D.h"
#include "Rn.h"

float *p2o3D,*p2s3D,*work3D1,*work3D2;
int stereo3D;
float *penpos3Ds;
float dist3D,sep3D;


int Init3DPlot(int stereo,float dist,float sep) {
	float alpha;

	p2o3D=allocM(3,3);
	p2s3D=allocM(3,3);
	work3D1=allocM(3,3);
	work3D2=allocM(3,3);
	penpos3Ds=allocV(2);
	if(!p2o3D||!p2s3D||!work3D1||!work3D2||!penpos3Ds) return 1;
	InitPlot();
	stereo3D=stereo;
	dist3D=dist;
	sep3D=sep;
	if(stereo) SetScales(-1,3,-1,1);
	else SetScales(-1,1,-1,1);
	setstdM(p2o3D,3,3,1);
	setstdM(p2s3D,3,3,1);
	if(stereo) {
		alpha=atan(sep/(2*dist));
		Rotate3D(29,2*alpha);
		setstdM(p2s3D,3,3,1);
		Rotate3D(28,alpha); }
	setstdV(penpos3Ds,2,0);
	return 0; }

Pt3Dptr Pt3DAlloc() {
	Pt3Dptr pt;

	pt=(Pt3Dptr) malloc(sizeof(struct Point3D));
	if(!pt) return NULL;
	pt->r[0]=pt->r[1]=pt->r[2]=0;
	pt->p[0]=pt->p[1]=pt->p[2]=0;
	pt->o[0]=pt->o[1]=pt->o[2]=0;
	pt->s[0]=pt->s[1]=pt->s[2]=0;
	pt->r2p=NULL;
	return pt; }

void Pt3DFree(Pt3Dptr pt) {
	if(pt) free(pt);
	return; }

void PlotPt3D(Pt3Dptr pt) {
	PlotPt(pt->o[0],pt->o[1]);
	if(stereo3D) {
		PlotPt(pt->s[0]+2,pt->s[1]);
		penpos3Ds[0]=pt->s[0]+2;
		penpos3Ds[1]=pt->s[1];
		PlotMove(pt->o[0],pt->o[1]); }
	return; }

void PlotLine3D(Pt3Dptr pt) {
	PlotLine(pt->o[0],pt->o[1]);
	if(stereo3D) {
		PlotMove(penpos3Ds[0],penpos3Ds[1]);
		PlotLine(pt->s[0]+2,pt->s[1]);
		penpos3Ds[0]=pt->s[0]+2;
		penpos3Ds[1]=pt->s[1];
		PlotMove(pt->o[0],pt->o[1]); }
	return; }

void PlotStr3D(Pt3Dptr pt,char *s) {
	PlotStr(pt->o[0],pt->o[1],s);
	if(stereo3D) {
		PlotStr(pt->s[0]+2,pt->s[1],s);
		penpos3Ds[0]=pt->s[0]+2;
		penpos3Ds[1]=pt->s[1];
		PlotMove(pt->o[0],pt->o[1]); }
	return; }

void Rotate3D(char c,float alpha) {
	float *m,*temp;

	if(!c) return;
	m=work3D1;
	setstdM(m,3,3,1);
	if(c==31) {m[4]=m[8]=cos(alpha);m[5]=-(m[7]=sin(alpha)); }
	else if(c==30) {m[4]=m[8]=cos(alpha);m[5]=-(m[7]=-sin(alpha)); }
	else if(c==28) {m[0]=m[8]=cos(alpha);m[2]=-(m[6]=sin(alpha)); }
	else if(c==29) {m[0]=m[8]=cos(alpha);m[2]=-(m[6]=-sin(alpha)); }
	else if(c=='.') {m[0]=m[4]=cos(alpha);m[1]=-(m[3]=sin(alpha)); }
	else if(c=='/') {m[0]=m[4]=cos(alpha);m[1]=-(m[3]=-sin(alpha)); }
	else return;
	dotMM(m,p2o3D,work3D2,3,3,3);
	temp=p2o3D;p2o3D=work3D2;work3D2=temp;
	if(stereo3D) {
		dotMM(m,p2s3D,work3D2,3,3,3);
		temp=p2s3D;p2s3D=work3D2;work3D2=temp; }
	return;	}

void UpdatePt3D(Pt3Dptr pt,int internal) {
	if(internal&&pt->r2p) dotVM(pt->r,pt->r2p,pt->p,3,3);
	dotMV(p2o3D,pt->p,pt->o,3,3);
	pt->o[0]*=dist3D/(dist3D-pt->o[2]);
	pt->o[1]*=dist3D/(dist3D-pt->o[2]);
	if(stereo3D) {
		dotMV(p2s3D,pt->p,pt->s,3,3);
		pt->s[0]*=dist3D/(dist3D-pt->s[2]);
		pt->s[1]*=dist3D/(dist3D-pt->s[2]); }
	return; }


Draw3Dptr Draw3DAlloc(int max,float *r2p) {
	Draw3Dptr dr;
	Pt3Dptr pt;
	int i;
	
	dr=(Draw3Dptr) malloc(sizeof(struct Drawing3D));
	if(!dr) return NULL;
	dr->max=max;
	dr->n=0;
	dr->pt=NULL;
	dr->col=NULL;
	dr->line=NULL;
	dr->r2p=r2p;
	dr->pt=(Pt3Dptr *) calloc(sizeof(Pt3Dptr),max);
	if(dr->pt) for(i=0;i<max;i++) dr->pt[i]=NULL;
	dr->col=(char *) calloc(sizeof(char),max);
	dr->line=(int *) calloc(sizeof(int),max);
	if(!dr->pt||!dr->col||!dr->line) {Draw3DFree(dr);return NULL;}
	for(i=0;i<max;i++) {
		pt=dr->pt[i]=Pt3DAlloc();
		if(!pt) {Draw3DFree(dr);return NULL;}
		pt->r[0]=pt->r[1]=pt->r[2]=0;
		pt->p[0]=pt->p[1]=pt->p[2]=0;
		pt->o[0]=pt->o[1]=pt->o[2]=0;
		pt->s[0]=pt->s[1]=pt->s[2]=0;
		pt->r2p=r2p;
		dr->col[i]='0';
		dr->line[i]=0; }
	return dr; }

void Draw3DFree(Draw3Dptr dr) {
	int i;

	if(!dr) return;
	if(dr->pt) for(i=0;i<dr->max;i++) Pt3DFree(dr->pt[i]);	
	if(dr->pt) free(dr->pt);
	if(dr->col) free(dr->col);
	if(dr->line) free(dr->line);
	free(dr);
	return; }

int DrawingPt3D(Draw3Dptr dr,float x,float y,float z,int rel,int line,int col) {
	int n;
	Pt3Dptr pt;

	if(dr->n==dr->max) return 1;
	n=dr->n++;
	pt=dr->pt[n];
	if(rel) {pt->r[0]=x;pt->r[1]=y;pt->r[2]=z;}
	else {pt->p[0]=x;pt->p[1]=y;pt->p[2]=z;}
	pt->r2p=dr->r2p;
	dr->col[n]=col;
	dr->line[n]=line;
	return; }

void PlotDraw3D(Draw3Dptr dr) {
	int i;

	for(i=0;i<dr->n;i++) {
		SetColor(dr->col[i]);
		if(dr->line[i]) PlotLine3D(dr->pt[i]);
		else PlotPt3D(dr->pt[i]); }
	return; }

void UpdateDraw3D(Draw3Dptr dr,int internal) {
	int i;
	
	for(i=0;i<dr->n;i++) UpdatePt3D(dr->pt[i],internal);
	return; }

float *Draw3DOut(Draw3Dptr dr,char c) {
	int i,j;
	float *a;
	
	a=allocM(dr->n,3);
	if(!a) return NULL;
	for(i=0;i<dr->n;i++)
		for(j=0;j<3;j++)
			if(c=='r') a[3*i+j]=dr->pt[i]->r[j];
			else if(c=='p') a[3*i+j]=dr->pt[i]->p[j];
			else if(c=='s') a[3*i+j]=dr->pt[i]->s[j];
			else a[3*i+j]=dr->pt[i]->o[j];
	return a; }

int Draw3D2ps(Draw3Dptr dr,int zdiv,char *name,int append) {
	FILE *fptr;
	char str[256];
	Pt3Dptr *pt;
	int nmfre,er,i,div,divcur;
	float zmin,zmax;
	
	er=GetName(&name,&nmfre);
	if(er) return er;
	if(!strcmp(name,"stdout")) fptr=stdout;
	else if(append==1) fptr=fopen(name,"a");
	else if(append==2) fptr=fopen(name,"w");
	else {
		fptr=fopen(name,"r");
		if(fptr) {
			fclose(fptr);
			fprintf(stderr,"File exists. Do you want to overwrite it? ");
			scanf("%s",str);
			if(str[0]!='y'&&str[0]!='Y') {
				if(nmfre) free(name);
				return 2;}}
		fptr=fopen(name,"w");}
	if(nmfre) free(name);
	if(!fptr) return 3;

	pt=dr->pt;
/*	zmin=zmax=pt[0]->o[2];
	for(i=0;i<dr->n;i++) {
		if(pt[i]->o[2]<zmin) zmin=pt[i]->o[2];
		else if(pt[i]->o[2]>zmax) zmax=pt[i]->o[2]; }*/
	zmin=-1;
	zmax=1;

	fprintf(fptr,"%!\n");
	fprintf(fptr,"%%%% Postscript file produced by Plot3D.c\n");
	fprintf(fptr,"/unit {250 mul} def\n");
	div=divcur=floor((pt[0]->o[2]-zmin)/(zmax-zmin)*zdiv);
	if(div<0) div=0;
	else if(div>=zdiv) div=zdiv-1;
	fprintf(fptr,"newpath\n");
	fprintf(fptr,"%f setlinewidth\n",0.5*(1.0-1.0*div/zdiv)-0.3);
	fprintf(fptr,"%f unit %f unit moveto\n",3+pt[0]->o[0],6+pt[0]->o[1]);
	for(i=1;i<dr->n;i++) {
		div=floor((zmax-pt[i]->o[2])/(zmax-zmin)*zdiv);
		if(div<0) div=0;
		else if(div>=zdiv) div=zdiv-1;
		if(div==divcur&&dr->line[i])
			fprintf(fptr,"%f unit %f unit lineto\n",3+pt[i]->o[0],6+pt[i]->o[1]);
		else {
			if(dr->line[i])
				fprintf(fptr,"%f unit %f unit lineto\n",3+pt[i]->o[0],6+pt[i]->o[1]);
			fprintf(fptr,"stroke\n");
			fprintf(fptr,"newpath\n");
			fprintf(fptr,"%f setlinewidth\n",(1.0-1.0*div/zdiv)-0.3);
			divcur=div;
			fprintf(fptr,"%f unit %f unit moveto\n",3+pt[i]->o[0],6+pt[i]->o[1]); }}
	fprintf(fptr,"stroke\n");
	fprintf(fptr,"showpage\n");
	if(fptr!=stdout) fclose(fptr);
	return er; }