/* Sample Input File: sf.in 
4000
0 0
10000 0 0 4
10 5
1
10
-5 3 -10
-5
1
-10
5
3 */
/*
Copyright (C) 2003  Mihai Pruna

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
http://www.gnu.org/licenses/gpl.html

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <conio.h>
#include <math.h>
float Ix,Iy,Ixy,vx,vy,gx,gy,denom,distance, height,x,y,vxx,vxy,vyx,vyy;
float q[50], sumad, A, diff [50], area[50], pprime[50],sx[50], sy[50], sa[50];
int n,i;
FILE *input;


void load() {

	input=fopen("sf.in","rt");


fscanf (input,"%f",&vx);



fscanf (input,"%f",&vxx);



fscanf (input, "%f",&vxy);


fscanf (input, "%f", &vy);

fscanf (input, "%f",&vyx);


fscanf (input, "%f",&vyy);


fscanf (input, "%d", &n);


	for (i=1;i<=n;i++) {


	fscanf (input," %f", &sx[i]);

	fscanf (input, " %f", &sy[i]);

	fscanf (input, " %f", &sa[i]);
	}



	fclose (input);
	};


void minput()
{

printf ("\n Input horiz. load " );

scanf ("%f",&vx);


printf ("\n Input x location of horiz. load ");

scanf ("%f",&vxx);


printf ("\n Input y location of horiz. load ");

scanf ("%f",&vxy);



printf ("\n Input vert. load ");

scanf ("%f", &vy);

printf ("\n Input x location of vert. load ");

scanf ("%f",&vyx);


printf ("\n Input y location of vert. load ");

scanf ("%f",&vyy);

printf ("\n Input number of stringers " );

scanf ("%d", &n);

printf ("\n Input stringers-clockwise");


	for (i=1;i<=n;i++) {

	printf ("\n Input x loc. for stringer %d ", i );

	scanf (" %f", &sx[i]);

	printf ("\n Input y loc. for stringer %d ",i);

	scanf (" %f", &sy[i]);

	printf ("\n Input area for stringer %d ", i );

	scanf (" %f", &sa[i]);
	}


};




void main()
{


char c;

printf ("\n Shear flows calculator for an idealized closed box structure");
printf ("\n using methods taught by Dr. Frank Radosta and Dr. Habib Eslami");
printf ("\n Copyright 2003 Mihai Pruna");
printf ("\n !!!This code has been tested for accuracy, however, not all particular cases   have been studied.");
printf ("\n The author assumes no responsibility for the accuracy of the results obtained  using this code!!!");


printf ("\n Select option: m-Manual Input; l-Load Input from file");

c=getch();
switch (c){

	case 'm':
		minput ();
	break;
	case 'l':
		load();
	break;
}





printf ("\n Now calculating centroid");

gx=0;
gy=0;
denom=0;


for (i=1;i<=n;i++)

{
	gx+=sx[i]*sa[i];
	denom+=sa[i];
	gy+=sy[i]*sa[i];
}
gx=gx/denom;
gy=gy/denom;

printf ("\nCentroid at %3.3f,%3.3f   Press any key to continue",gx,gy);


getch();





printf ("\n Recentering around centroid...");

for (i=1;i<=n;i++){
	sx[i]=sx[i]-gx;
	sy[i]=sy[i]-gy;
	printf ("\n Stringer %d located at %3.3f, %3.3f", i,sx[i], sy[i]);
}



vxx=vxx-gx;
vyx=vyx-gx;
vxy=vxy-gy;
vyy=vyy-gy;

printf ("\nVertical load of %3.3f at %3.3f,%3.3f, Horizontal load of %3.3f at %3.3f,%3.3f",vy,vyx,vyy,vx,vxx,vxy );




printf ("\n Now calculating moments of inertia");

Ix=0;
Iy=0;
Ixy=0;
for (i=1;i<=n;i++){

Ix+=sy[i]*sy[i]*sa[i];
Iy+=sx[i]*sx[i]*sa[i];
Ixy+=sx[i]*sy[i]*sa[i];
}

printf ("\n Ix is %3.3f, Iy is %3.3f, Ixy is %3.3f  Press any key", Ix, Iy, Ixy);
getch();

printf ("\n Now calculating P primes");

for (i=1;i<=n;i++){

pprime[i]=0.000;
	pprime[i]=(-(1/(Ix*Iy-Ixy*Ixy))*((Iy*vy-Ixy*vx)*sy[i]+(Ix*vx-Ixy*vy)*sx[i])*sa[i]);

printf ("\n Pprime (%d) is %3.3f ",i, pprime[i] );

}





printf ("\n Calculating areas");


for (i=2;i<=n;i++) {

	distance=sqrt(sx[i]*sx[i]+sy[i]*sy[i]);
x=(sy[i-1]+(sx[i-1]*sx[i]/sy[i]))/(sy[i]/sx[i]+sx[i]/sy[i]);
y=x*sy[i]/sx[i];
height=sqrt((sx[i-1]-x)*(sx[i-1]-x)+(sy[i-1]-y)*(sy[i-1]-y));
area[i-1]=distance*height/2;
printf ("\n Area %d is %3.3f",i-1,area[i-1]);
}


	distance=sqrt(sx[1]*sx[1]+sy[1]*sy[1]);
x=(sy[n]+(sx[n]*sx[1]/sy[1]))/(sy[1]/sx[1]+sx[1]/sy[1]);
y=x*sy[1]/sx[1];
height=sqrt((sx[n]-x)*(sx[n]-x)+(sy[n]-y)*(sy[n]-y));
area[n]=distance*height/2;
printf ("\n Area %d is %3.3f",n,area[n]);

printf ("\n Press any key to continue");
getch();

diff[1]=pprime[2];

printf ("\n q2-q1=%3.3f",diff[1]);
for (i=2;i<=n-1;i++)
	{
		diff[i]=diff[i-1]+pprime[i+1];

		printf ("\n q%d-q1=%3.3f",i+1,diff[i]);

}

A=0;
printf ("\n Calculating shear flows. Clockwise Direction is assumed positive");
printf("\n Press any key");
getch ();

for (i=1;i<=n;i++){

	A+=area[i];}

	sumad=0;

printf ("\n Total area is %3.3f ",A);

	for (i=1;i<=n-1;i++) {
		sumad+=2*diff[i]*area[i+1];}

q[1]=0.000;
		q[1]=(vx*vxy+vy*vyx-sumad)/(2*A);

		printf("\n q1 is %3.3f ",q[1]);

for (i=2;i<=n;i++){
q[i]=0.000;
q[i]=q[1]+diff[i-1];
printf ("\n q%d is %3.3f ",i,q[i]);
}

printf ("\n Press any key to end program");
getch();



}