[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Isoclines



> If I understand correctly, you want to draw level curves of f(x,y).  One
> trick that's worked for me is to compute J(grad f), the gradient of f
> rotated a quarter turn, then to draw solution curves of the resulting
> planar ODE. It's not elegant; 
Thank you for your solution, but if that's not elegant, then what is my 
solution included below? Ugly?

#include "epix.h"
#include <cstdlib>
using namespace ePiX;
double ISOKLINC=1; // must unfortunatly use global variable

double f(double x){
	if (ISOKLINC!=0.0)
		return -x/ISOKLINC;
	else
		return 10000;
}

// This function uses the global variable ISOKLINC and f uses it also.
void isoklin(double f(double x), double minx,double maxx,double steps,double 
k)
{
	ISOKLINC=k;
	double h=(maxx-minx)/steps;
	P p;
	P kp=P(1,ISOKLINC);
	kp=0.8/norm(kp) * kp;
	for (double t=minx;t<=maxx;t=t+h)
	{
		p=P(t,f(t));// the point
		line(p-kp,p+kp);
	}
	char c[20];
	int a=ISOKLINC;
	if (a==ISOKLINC)
		sprintf(c,"%i",a);
	else
		sprintf(c,"%8.2f",ISOKLINC);
	label(P(3+maxx+0.2,f(maxx+0.2)),"$C="+std::string(c)+"$");
	plot(f,minx,maxx,30);
}

int main(){
	unitlength("1cm");
	picture(5,5);
	double x1=-10;
	double x2=10;
	double y1=-10;
	double y2=10;
	double N=(x2-x1)/2.0;
	bounding_box(P(x1,y1),P(x2,y2));

	double x,y,k;
	begin();
		dot_size(1);
		pen(0.2);
		isoklin(f,x1,x2,N,2);
		isoklin(f,x1,x2,N,-2);
		isoklin(f,x1,x2,N,-1);
		isoklin(f,x1,x2,N,1);
		isoklin(f,x1,x2,N,-8);
		isoklin(f,-2.4,2.4,N,0.-0.21);
		isoklin(f,-3.6,3.6,N,0.31);
		pen(1);
		h_axis(1);
		v_axis(1);
	end();
}