#! /usr/bin/env python

import os, string, sys, cp, commands

mandelSource = r"""

                /********************************/
                /*                              */
                /*     The Mandelbrot Set       */
                /*                              */
                /*      z = z*z + c             */
                /*                              */
                /********************************/
                
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

#define LARGE_RADIUS 2.0
                                        /* Maximum number of iterations */
                                        /* before declaring a point in  */
                                        /* the Mandelbrot set           */
#define MAX_ITERATIONS 100
#define WIDTH 100

int main(int argc, char *argv[]){
	
	double z_real, z_img, z_magnitude;
	double c_real, c_img, crmin, crmax, cimin, cimax, dcr, dci, z_current_real;
	int i, j;/*, NX, NY;*/
	int counter;
	int val;
	int nx, ny, cycles;
	double scalefactor = (double)(255) / (double)MAX_ITERATIONS;
	FILE *fout;
	
	if (argc < 2)
	{
		printf("Usage: rmin rmax imin imax gridwidth cycles\n");
		return 0;
	}
	
/*	printf("is 1st arg same as def?: %i\n", strcmp(argv[1],"def"));*/
	if (strcmp(argv[1],"def") == 0)
	{
		/* define grid spanning the complex plane for the C variable */
		crmin= -1.7;
		crmax = 0.8;
		cimin = -1.0;
		cimax = 1.0;
		nx = WIDTH;
		ny = WIDTH;
		cycles = MAX_ITERATIONS;
		printf("default values: \nrmin= %f rmax= %f \nimin= %f imax= %f \nwidth= %i cycles= %i\n", 
		crmin, crmax, cimin, cimax, WIDTH, cycles);
	}
	else
	{
		/* define grid spanning the complex plane for the C variable */
		
		crmin = atof(argv[1]);
		crmax = atof(argv[2]);
		cimin = atof(argv[3]);
		cimax = atof(argv[4]);
		nx = atoi(argv[5]);
		ny = atoi(argv[5]);
		cycles = atoi(argv[6]);
		/*note: filename is in argv[6]*/
		printf("input values:\nrmin= %f rmax= %f \nimin= %f imax= %f \nwidth= %i cycles= %i\n", 
		crmin, crmax, cimin, cimax, nx, cycles);	
	}

	dcr = (crmax - crmin)/(nx-1);
	dci = (cimax - cimin)/(ny-1);
	
	/*open the file*/
	/*fout = fopen("outfile", "w");*/
	fout = fopen(argv[7], "wb");
	
	for ( i=0 ; i<nx; i++ )
	{
		for ( j=0 ; j<ny; j++ )
		{
			c_real = crmin + i*dcr;
			c_img = cimin + j*dci;
			/* z starts at origin */
			z_real = 0.0;
			z_img  = 0.0;
			/* set counter to zero */
			counter = 0;
			/* iterate map for MAX_ITERATIONS times or ...*/
			while ( counter < cycles ) 
			{
				z_current_real = z_real;
				z_real = z_real*z_real - z_img*z_img + c_real;
				z_img  = 2.0*z_current_real*z_img + c_img;
				counter = counter + 1;
				/* ... until magnitude of z */
				/* is larger than some      */
				/* large radius             */
				z_magnitude = z_real*z_real+z_img*z_img;
				if ( z_magnitude > LARGE_RADIUS ) break;
			}
			
			/* write iteration number for all scanned C values for color picture  */
			
			val = (int)floor(counter * scalefactor);
/*			fprintf ( fout, "%i ", val);*/
/*			fwrite( val, 8, 1, fout);*/
			fputc(val, fout);

		}
	
	}
	fclose( fout );

	return 0;
}
"""


def main():
	"""
mandel.py

Calls mandel.c to create an array of numbers generated 
using the Mandelbrot formula. Attempts to compile mandel.c 
if mandel isn't found on PATH. Writes an 8-bit binary file which
can be viewed in NIH Image or Photoshop, for example.

c		comma-delimited list of complex plane 
		delimiters: rmin,rmax,imin,imax
w		width of grid
n		number of cycles
v 		verbosity
nc		new coordinates - give pixel coords corresponding to some
		area within the plot of width w, bounded by coordinates c;
		will return new values for c. 
		new coordinates 
out		name of outfile

NOTE: To import into NIH Image, choose File --> Import
and select Custom, 8-bit, and set to appropriate grid width.

Rotate -90 degrees to maintain proper orientation!  
"""
	
	
	optlist = [('c','-1.7,0.8,-1.0,1.0'),
	('w',1000),
	('n',100),
	('v',0),
	('h',),
	('out','out.mandel'),
	('nc','') ]
	
	dict = cp.commandparser(options=optlist, usage='-h for usage')
	if dict.status('h'):		
		print main.__doc__
		sys.exit(0)
		
	foundit = findSource('mandel', 0)
	
	if not foundit: #compile executable from code above
		f = open('mandel.c', 'w')
		f.write(mandelSource)
		f.close()
		
		# try to compile
		print 'mandel executable not found: compiling from mandel.c'
		tup = commands.getstatusoutput( 'cc -o mandel mandel.c' )
		if tup[0] == 0: # exit status 0 = ok
			print tup[1]
		else:
			sys.exit('compilation failed:\n', cmd)
	
	coordinates = dict.value('c').split(',')
	
	try:
		rmin, rmax, imin, imax = float(coordinates[0]), \
		float(coordinates[1]), \
		float(coordinates[2]), \
		float(coordinates[3])
	except IndexError: 
		sys.exit('-c must be a comma-delimited list of 4 values')
	
	gridSize = dict.value('w')
	cycles = dict.value('n')
	v = dict.value('v')
	outfilename = dict.value('out')
	newCoords = dict.value('nc')
	
	if newCoords != '':
		try:
			newCoords = newCoords.split(',')
			rmin_pix, rmax_pix, imin_pix, imax_pix = float(newCoords[0]), \
			float(newCoords[1]), \
			float(newCoords[2]), \
			float(newCoords[3])
		except IndexError: 
			sys.exit('-c must be a comma-delimited list of 4 values')
			
		rmin_new = rmin + float(rmin_pix)*(rmax-rmin)/gridSize
		rmax_new = rmin + float(rmax_pix)*(rmax-rmin)/gridSize
		imin_new = imin + float(imin_pix)*(imax-imin)/gridSize
		immax_new = imin + float(imax_pix)*(imax-imin)/gridSize
		
		sys.exit('New coords: %f,%f,%f,%f' % (rmin_new, rmax_new, imin_new, immax_new) )
		
	
	print 'Writing', outfilename
	
	cmd = 'mandel %f %f %f %f %i %i %s' % (rmin, rmax, imin, imax, gridSize, cycles, outfilename)
	
	tup = commands.getstatusoutput( cmd )
	if tup[0] == 0: # exit status 0 = ok
		print tup[1]
	else:
		print 'failed:', cmd
		
def findSource(filename, v):
	pathstr = commands.getstatusoutput( "echo $PATH" )[1]
	pathlist = string.split(pathstr, ':')
	outpath = ''
	for location in pathlist:
		if location[-1] != "/":
			location = location + "/"
		path = location + filename
		if v: print location
		if os.access(path, os.R_OK):
			outpath = path
			break
	found = 1
	if outpath == '':
		found = 0
	return found
	
if __name__ == '__main__':
	main()
	#testGapNumbering()