Subject: Re: CCL:calculate the angle

From: [email-masked] (colvin michael e)
Date: Mon, 26 Feb 1996 14:29:03 -0800

Recently [email-masked] (Pornthep Sompornpisut) asked:

>I will appreciate if anyone know of any programs that can measure the angle
>between the planes? 

You are probably looking for something that has a nice GUI for selecting
the atoms, but if you just want a program to calculate the angles, 
here's a tiny C code I wrote sometime back to do this.  It just reads sets
of six cartesian points, three each per plane, and prints out the angle.
Nota bene: To the best of my knowledge it works, but of course I offer
no guarantees to its accuracy, reliability, suitability, etc. etc..., so
definitely check this against a hand-calculated result.

                                               --Mike Colvin
                                                 [email-masked]

=============================================================
/* Routine to calculate the angle between two planes */

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

typedef struct {
	double x;
	double y;
	double z;
	} point;
	
int main(void);
point cross(point, point);
double dot(point, point);

main()
{
    point a1,a2,a3,b1,b2,b3;
    point va1,va2,vb1,vb2, crossa, crossb;
    double dotab, angle;
	
    while(1)
    {
        if (scanf("%lf %lf %lf",&a1.x,&a1.y,&a1.z)==EOF)
            exit(0);
        scanf("%lf %lf %lf",&a2.x,&a2.y,&a2.z);
        scanf("%lf %lf %lf",&a3.x,&a3.y,&a3.z);
        scanf("%lf %lf %lf",&b1.x,&b1.y,&b1.z);
        scanf("%lf %lf %lf",&b2.x,&b2.y,&b2.z);
        scanf("%lf %lf %lf",&b3.x,&b3.y,&b3.z);
		
        /* calculate direction vectors */
        va1.x=a2.x-a1.x; va1.y=a2.y-a1.y; va1.z=a2.z-a1.z;
        va2.x=a3.x-a1.x; va2.y=a3.y-a1.y; va2.z=a3.z-a1.z;
        vb1.x=b2.x-b1.x; vb1.y=b2.y-b1.y; vb1.z=b2.z-b1.z;
        vb2.x=b3.x-b1.x; vb2.y=b3.y-b1.y; vb2.z=b3.z-b1.z;
		
        /* calculate cross products */
        crossa=cross(va1,va2);
        crossb=cross(vb1,vb2);
		
        dotab=dot(crossa,crossb);
        dotab/=sqrt(dot(crossa,crossa));
        dotab/=sqrt(dot(crossb,crossb));
		
        angle=acos(dotab)*180/3.1415926535;
		
        printf(" Angle = %lf\n",angle);
    }
}

point cross(point a, point b)
{
	point value;
	
	value.x=a.y*b.z-a.z*b.y;
	value.y=a.z*b.x-a.x*b.z;
	value.z=a.x*b.y-a.y*b.x;
	
	return value;
}

double dot(point a, point b)
{
	return a.x*b.x+a.y*b.y+a.z*b.z;
}
=============================================================