Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

More Problems : 2d array multiplication 1

Status
Not open for further replies.

aaadetos

Programmer
Nov 21, 2004
54
US
Can someone please tell me why vargm() is crashing in main()?
Code:
// Q3.cpp : Defines the entry point for the console application.
//

#include <math>
#include "stdafx.h"
#include <iostream>
using namespace std;

double t, p;	//azimuth and dip angles respectively
double lambda1 = 2.2;
double lambda2 = 1;//hz = 0, so this will have no influence for the 0 degrees dip case
double a11,a12,a21,a22,d11,d13,d31,d33;
double lagh,hy,hx,Y1,Y2,C0,Y,y_1,y_2,var1,var2;
double ax = 55;// range of variogram in E-W direction: fixed
double ay = 25;//range of variograms in N-S directions: fixed 

double vargm();

const double pi = 3.141592653589793238462643383279502884197;
const int rows = 3;
const int cols = 3;

double lag,f_lag;

//---------------------------------------------------------
//Initializing the arrays to zero
//---------------------------------------------------------
	//array that holds the dip angle, phi
	double dip_array [3][3] = {0};	

	//array that holds the azimuth angle, theta
	double azim_array [3][3] = {0};	

	//array that holds the anisotropy ratios
	double lambda_array [3][3] = {0};	

	//array that holds the original lag values
	double r_array [3] = {0};

	//array that holds final lag values
//	double fr_array [3][3] = {0};	

	//array that holds final lag values
	double fr_array2 [3] = {0};	



//---------------------------------------------------------
//Code to convert the angles to radians, acceptable in C++
//---------------------------------------------------------
//Dip Angle,p
double a(double p)
{
	d11 = cos((p*pi)/180);
	return d11;
}

double b(double p)
{
	d13 = -1*sin((p*pi)/180);
	return d13;
}

double c(double p)
{
	d31 = sin((p*pi)/180);
	return d31;
}

double d(double p)
{
	d33 = cos((p*pi)/180);
	return d33;
}


//Azimuth Angles, t
double e(double t)
{
	a11 = sin((t*pi)/180);
	return a11;
}

double f(double t)
{
	a12 = cos((t*pi)/180);
	return a12;
}

double g(double t)
{
	a21 = -1*cos((t*pi)/180);
	return a21;
}

double h(double t)
{
	a22 = sin((t*pi)/180);
	return a22;
}


void printarray (double arg[][3]) 
{
  for (int i=0; i<rows; i++)//for each row
  {
	  for (int j=0; j<cols; j++)	//for each column value
		  cout << arg[i][j] << " ";
	  cout << "\n";
  }
}

void printarray2(double arg[]) 
{
  for (int i=0; i<rows; i++)//for each row
  {
	  //for (int j=0; j<3; j++)	//for each column value
		  cout << arg[i] << " ";
	  cout << "\n";
  }
}

int main(int argc, char* argv[])
{
	
	cout<<"ANISOTROPY MODELLING.\n";
	cout<<"\nPlease enter a real number between 0 and 360 for the azimuth angle, theta, then press the 'Enter' key.\n\n";
	cin>>t;

	cout<<"\nPlease enter a real number between 0 and 360 for the dip angle, phi, then press the 'Enter' key.\n\n";
	cin>>p;

	a(p);
	b(p);
	c(p);
	d(p);
	e(t);
	f(t);
	g(t);
	h(t);

	//-----------------------------------------------------------
	//Loading the matrices
	//-----------------------------------------------------------
	//array that holds the dip angle, phi
	double dip_array [3][3] = {{d11,0,d13},{0,1,0},{d31,0,d33}};	

	//array that holds the azimuth angle, theta
	double azim_array [3][3] = {{a11,a12,0},{a21,a22,0},{0,0,1}};	

	//array that holds the anisotropy ratios
	double lambda_array [3][3] = {{lambda1,0,0},{0,1,0},{0,0,lambda2}};	

	//array that holds the original lag values
	cout<<"\nPlease enter the lag distance in the E-W (x-direction) and then press the 'Enter' key.\n\n";
	cin>>hx;
	cout<<"\nhx = "<<hx<<endl;//"  "<<"ax = " <<ax<<endl;

	cout<<"\nPlease enter the lag distance in the N-S (y-direction) and then press the 'Enter' key.\n\n";
	cin>>hy;
//	cout<<"\nPlease enter the range in the N-S (y-direction) and then press the 'Enter' key.\n\n";
//	cin>>ay;
	cout<<"\nhy = "<<hy<<endl;//"  "<<"ay = "<<ay<<endl;

	double r_array [3] = {hx,hy,0};


	//----------------------------------------------------------------

	//code to check that the arrays are properly initialized in the matrices for any theta,phi
	cout<<"\nThe dip array is:\n";
	printarray(dip_array);

	cout<<"\nThe azimuth array is:\n";
	printarray(azim_array);

	cout<<"\nThe anisotropy array is:\n";
	printarray(lambda_array);

	cout<<"\nThe original lag array is:\n";
	printarray2(r_array);

	//----------------------------------------------------------------

	//Matrix multiplication code
	for(int i=0;i<rows;i++)	//for each row
		{
			for(int j=0;j<cols;j++)	//for each column
			{
				for(int k=0;k<cols;k++)	//for each column
				{
					for(int l=0;l<cols;l++)	//for each column
					{
						fr_array2 [i] += lambda_array[i][l]*dip_array[l][k]*azim_array[k][j]*r_array[j];
					}
				}
			}
		}
	
	//Print the final lag array	
	cout<<"\nThe final lag array is:\n";
	printarray2(fr_array2);

	//code for determining the value of the lag
	for(i=0;i<rows;i++)
	{
		lag += pow(fr_array2[i],2);
		f_lag = pow(lag,0.5);
	}


	//Print the final lag value
	cout<<"\nThe final lag value - given the specified azimuth and dip angles - is:\n";
	cout<<"lag = "<<f_lag;//Prgm crashes beyond this point!! HEEELP!!!

	vargm();

	//Print the final variogram value
	cout<<"\n\nThe value of the variogram - given the specified azimuth and dip angles - is:\n";
	cout<<"Y(h) = "<<Y<<endl;


	return 0;
}

//code for spherical variograms
double vargm()
{
	if (t == 0 && p == 0)
	{
		lagh = sqrt(pow(hx,2)+pow(hy,2));
		while (lagh<ay) //lagh/ay<1
		{
			var1 = 1.5*sqrt(pow(lagh/ax,2));
			var2 = 0.5*pow(pow(lagh/ax,2),3);
			y_1 = var1 - var2;
		}
		while (lagh/ay>1)
		{
			y_1 = 1;
		}
		while (hx/ax<1)
		{
			//	h = sqrt(pow(hx,2)+pow(hy,2));
			var1 = 1.5*sqrt(pow(hx/ax,2));
			var2 = 0.5*pow(pow(hx/ax,2),3);
			y_2 = var1 - var2;
		}
		while (hx/ax>1)
		{
			y_2 = 1;
		}

	}

	if (t != 0 || p != 0)
	{
		lagh = sqrt(pow(hx,2)+pow(hy,2));
		while (lagh<ay)
		{
			hx = fr_array2[0];
			hy = fr_array2[1];
			var1 = 1.5*sqrt(pow((lagh/ax),2));
			var2 = 0.5*pow(pow((lagh/ax),2),3);
			y_1 = var1 - var2;
		}
		while (lagh>ay)
		{
			y_1 = 1;
		}
		while (hx/ax<1)
		{
			hx = fr_array2[0];
			hy = fr_array2[1];
			//	h = sqrt(pow(hx,2)+pow(hy,2));
			var1 = 1.5*sqrt(pow(hx/ax,2));
			var2 = 0.5*pow(sqrt(pow(hx/ax,2)),3);
			y_2 = var1 - var2;
		}
		while (hx/ax>1)
		{
			y_2 =1;
		}
	}
	
	C0 = 0.125 + 0.075;
	Y1 = 0.075*y_1;
	Y2 = 0.15*y_2;

	Y = C0 + Y1 + Y2;
	return Y;
}

 
It's a rather strange code in vargm(), for example:
Code:
while (lagh/ay>1)
{
   y_1 = 1;
}
// or
while (lagh<ay) //lagh/ay<1
{
   var1 = 1.5*sqrt(pow(lagh/ax,2));
   var2 = 0.5*pow(pow(lagh/ax,2),3);
   y_1 = var1 - var2;
}
Loop condition variables are intact in loop bodies. So these loops counts are zero or forever.

It's amusing to see this old good Fortran IV style (with unforgettable COMMONs) in C++...
 
Thanks, ArkM. I already replaced the while with if...else loops. Here's the new code:
Code:
double vargm()
{
	double h_x,h_y;
	lagh = sqrt(pow(hx,2)+pow(hy,2));
	if (t == 0 && p == 0)
	{
		if (lagh<ay) //lagh/ay<1
		{
			var1 = 1.5*sqrt(pow(lagh/ax,2));
			var2 = 0.5*pow(pow(lagh/ax,2),3);
			y_1 = var1 - var2;
		}
		else
		{
			y_1 = 1;
		}
		if (hx/ax<1)
		{
			var1 = 1.5*sqrt(pow(hx/ax,2));
			var2 = 0.5*pow(pow(hx/ax,2),3);
			y_2 = var1 - var2;
		}
		else
		{
			y_2 = 1;
		}

	}

	if (t != 0 || p != 0)
	{
		if(lagh<ay)
		{
			h_x = fr_array2[0];
			h_y = fr_array2[1];
			var1 = 1.5*sqrt(pow((lagh/ax),2));
			var2 = 0.5*pow(pow((lagh/ax),2),3);
			y_1 = var1 - var2;
		}
		else
		{
			y_1 = 1;
		}
		if (hx/ax<1)
		{
			h_x = fr_array2[0];
			h_y = fr_array2[1];
			//	h = sqrt(pow(hx,2)+pow(hy,2));
			var1 = 1.5*sqrt(pow(hx/ax,2));
			var2 = 0.5*pow(sqrt(pow(hx/ax,2)),3);
			y_2 = var1 - var2;
		}
		else
		{
			y_2 =1;
		}
	}
	
	C0 = 0.125 + 0.075;
	Y1 = 0.075*y_1;
	Y2 = 0.15*y_2;

	Y = C0 + Y1 + Y2;
	return Y;
}

 
Well, now it works ((in formal;).
Some little advice(s):
1. Don't use pow(x,2), better x*x (with possible intermediate vars as x = y/z; s = x * x;
2. I think it's more clear:
Code:
const double ppi  = p*pi/180;
const double pcos = cos(ppi);
const double psin = sin(ppi);
... // and so on, then
double array[3][3] = {{pcos,psin,....
instead of mysterious and cumbersome a(p), b(p) etc...
3. Of course, make a matrix multiplication function instead of implicit loops in main().
4. Stop global variables pollution. For example, lagh var used only in vargm(). It's a typical local variable.
5. You have in vargm:
Code:
    if (t == 0 && p == 0)
    { ....
    }

    if (t != 0 || p != 0)
    { ...
    }
Be more attentive: it's exactly:
Code:
if (t == 0 && p == 0)
   ...
else
   ...
It's much more clearer (and faster;)...
6. No comments:
Code:
// humorous
pow(pow(hx/ax,2),3);
// is
pow(hx/ax,6)
7. Use << endl instead of << "...\n" in user interface.
....
Good luck!
 
Thanks a lot...6. is real humour...tells u how "cautious" im thinking...lol!
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top