From 5b0492681b1ab58e15b6dda5e6196c93f4655dcd Mon Sep 17 00:00:00 2001
From: Vince Weaver <vincent.weaver@maine.edu>
Date: Wed, 1 Jun 2016 13:17:57 -0400
Subject: [PATCH] ksp: c: more work trying to get drag to work

---
 ksp/c/Makefile                    | 18 +++++++++--
 ksp/c/{atan_test.c => ksp_atan.c} |  0
 ksp/c/ksp_drag.c                  | 53 +++++++++++++++++++++++++++++++
 3 files changed, 68 insertions(+), 3 deletions(-)
 rename ksp/c/{atan_test.c => ksp_atan.c} (100%)
 create mode 100644 ksp/c/ksp_drag.c

diff --git a/ksp/c/Makefile b/ksp/c/Makefile
index 9a70f74a..21e64f2c 100644
--- a/ksp/c/Makefile
+++ b/ksp/c/Makefile
@@ -1,8 +1,8 @@
 CC = gcc
 CFLAGS = -Wall -O2
-LFLAGS =
+LFLAGS = -lm
 
-all:	ksp_acpx ksp_launch ksp_loading ksp_vab
+all:	ksp_acpx ksp_atan ksp_drag ksp_launch ksp_loading ksp_vab
 
 ksp_acpx:	ksp_acpx.o
 	$(CC) $(LFLAGS) -o ksp_acpx ksp_acpx.o
@@ -10,6 +10,18 @@ ksp_acpx:	ksp_acpx.o
 ksp_acpx.o:	ksp_acpx.c
 	$(CC) $(CFLAGS) -c ksp_acpx.c
 
+ksp_atan:	ksp_atan.o
+	$(CC) $(LFLAGS) -o ksp_atan ksp_atan.o
+
+ksp_atan.o:	ksp_atan.c
+	$(CC) $(CFLAGS) -c ksp_atan.c
+
+ksp_drag:	ksp_drag.o
+	$(CC) $(LFLAGS) -o ksp_drag ksp_drag.o
+
+ksp_drag.o:	ksp_drag.c
+	$(CC) $(CFLAGS) -c ksp_drag.c
+
 ksp_launch:	ksp_launch.o
 	$(CC) $(LFLAGS) -o ksp_launch ksp_launch.o -lm
 
@@ -30,5 +42,5 @@ ksp_vab.o:	ksp_vab.c
 
 
 clean:	
-	rm -f *~ *.o ksp_launch ksp_loading ksp_vab ksp_acpx
+	rm -f *~ *.o ksp_launch ksp_loading ksp_vab ksp_acpx ksp_drag ksp_atan
 
diff --git a/ksp/c/atan_test.c b/ksp/c/ksp_atan.c
similarity index 100%
rename from ksp/c/atan_test.c
rename to ksp/c/ksp_atan.c
diff --git a/ksp/c/ksp_drag.c b/ksp/c/ksp_drag.c
new file mode 100644
index 00000000..f0cf1f38
--- /dev/null
+++ b/ksp/c/ksp_drag.c
@@ -0,0 +1,53 @@
+#include <stdio.h>
+#include <math.h>
+
+int main(int argc, char **argv) {
+
+	double v0=0;
+	double g=9.8;
+	double mass=700;
+	double coeff=0.5;
+	double area=3.0*1000.0;
+	double density=0.1;
+	double time=0.0;
+	double a=0.0;
+	double altitude=20000.0;
+	double deltat=1.0;
+	double pressure,pressure0=101325.0;
+	double temperature=273.0;
+	double v_dir;
+	double ga,da;
+	double terminal_velocity;
+
+#define KERBIN_RADIUS	600000.0
+
+	while(1) {
+
+		terminal_velocity=sqrt( (2*mass*g)/(density*area*coeff));
+
+		if (v0>0) v_dir=-1; else v_dir=1;
+		ga=-g;
+		da=v_dir*(1.0/(2.0*mass))*coeff*density*area*v0*v0;
+
+		a=ga+da;
+
+//		v0=v0+(a*deltat);
+		v0=-terminal_velocity;
+		altitude=altitude+(v0*deltat);
+		time+=deltat;
+
+		pressure=pressure0*exp(-(altitude)/5600);
+                density=pressure/(287*temperature);
+
+		g=9.8/(
+                        ((altitude+KERBIN_RADIUS)/KERBIN_RADIUS)*
+                        ((altitude+KERBIN_RADIUS)/KERBIN_RADIUS));
+
+
+		printf("t=%lf ga=%lf da=%lf a=%lf v=%lf y=%lf tv=%lf\n",time,ga,da,a,v0,altitude,
+			terminal_velocity);
+		if ((altitude<0) || (altitude>21000)) break;
+	}
+
+	return 0;
+}