digitalredline's posterous

Gunnerkrigg Court, python style

So I'm a big fan of the web comic Gunnerkrigg Court, and I'm also a big fan of being able to read my comics like a book. Since I'm too cheap and lazy to buy the book, I decided to make one of my own (ish). With 734 comics (or something like that) currently online, downloading them one by one would have taken a LONG time. Probably much longer than it took me to cook up this script. After reading Toby Donaldson's Python Second Edition Visual Quickstart Guide, which I received from a friend, I more fully understood the power of Python scripts and desired to make one that did something useful. I also really, really want a iPad to use as a PDF and ebook reader, so this makes my comics way more iPad (and any portable media device) friendly. This program downloads all the files to a certain folder (The user specifies the path) and stops downloading when all of the files are downloaded. This program also runs as an updater that will count the number of comics and start downloading again from where the user's last comic is. That way, the updater can run every week and pick up three (hopefully) new comics. Please note that this script only works because the creator of the comic saves them as '00000000.jpg', so I just iterate through the numbers at the end of the eight digit string. I then used iCombiner to merge the jpg images into pdf documents (book one and book two). Note, these comics are property of the author and should not be used for commercial purposes.

iimport urllib

import os

# This code downloads and saves all of the comics available in the archives on www.gunnerkriggcourt.com to a specific directory '/file' of the users preference.

# By Mike McD

# Created: 6/14/10

# Last Updated: 6/15/10

# digitalredline.posterous.com

comicCounter=len(os.listdir('/file'))+1 # reads the number of files in the folder to start downloading at the next comic

errorCount=0

def download_comic(url,comicName):

"""

download a comic in the form of

url = http://www.example.com

comicName = '00000000.jpg'

"""

image=urllib.URLopener()

image.retrieve(url,comicName) # download comicName at URL

while comicCounter <= 1000: # not the most elegant solution

os.chdir('/file') # set where files download to

try:

if comicCounter < 10: # needed to break into 10^n segments because comic names are a set of zeros followed by a number

comicNumber=str('0000000'+str(comicCounter)) # string containing the eight digit comic number

comicName=str(comicNumber+".jpg") # string containing the file name

url=str("http://www.gunnerkrigg.com//comics/"+comicName) # creates the URL for the comic

comicCounter+=1 # increments the comic counter to go to the next comic, must be before the download in case the download raises an exception

download_comic(url,comicName) # uses the function defined above to download the comic

print url

if 10 <= comicCounter < 100:

comicNumber=str('000000'+str(comicCounter))

comicName=str(comicNumber+".jpg")

url=str("http://www.gunnerkrigg.com//comics/"+comicName)

comicCounter+=1

download_comic(url,comicName)

print url

if 100 <= comicCounter < 1000:

comicNumber=str('00000'+str(comicCounter))

comicName=str(comicNumber+".jpg")

url=str("http://www.gunnerkrigg.com//comics/"+comicName)

comicCounter+=1

download_comic(url,comicName)

print url

else: # quit the program if any number outside this range shows up

quit

except IOError: # urllib raises an IOError for a 404 error, when the comic doesn't exist

errorCount+=1 # add one to the error count

if errorCount>3: # if more than three errors occur during downloading, quit the program

break

else:

print str("comic"+ ' ' + str(comicCounter) + ' ' + "does not exist") # otherwise say that the certain comic number doesn't exist

print "all comics are up to date" # prints if all comics are downloaded

Click here to download:

GunnerKriggDL.py (2 KB)

Posted June 15, 2010 @ 12:54

// 0 Comments

Soar like an EAGLE

My newest toy, EAGLE, is making circuit diagrams and PBC design a breeze. Since I'm cheap and therefore using the "light" version, it has a few limitations (two sided, 3.2"x4"). Other than those, which I don't anticipate going over anytime soon. It's a lot of fun, but there is a rather steep learning curve associated with it. I've spent about an hour and a half on this guy (a 1.8V to 12V booster), getting the schematic laid up and then wiring it (since the auto-router function doesn't work in the free version).

Schematic:

Wiring Diagram:

Posted May 22, 2010 @ 9:59

// 0 Comments

Science Research Project: the halfway point.

/* Test program for tach lighting system  * By Mike McD  * Created on 5/8/2010.  Last Updated on 5/8/10.  * digitalredline.posterous.com  */  int LEDs[] = {0,22,23,2,3,4,5,6,7,8,9,10,11,12,13,14,15}; //Array of all LED pins int RPMs[] = {0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000}; //Array of all RPM values int RPMinput = 0; //Input for the pot int RPMvalue; //RPM value variable int RPM; //"True" rpm int n; //Index variable int remainder; //Remainder variable int PWMvalue; //PWM value variable  void setup() {   Serial.begin(9600);   pinMode(LEDs[1],OUTPUT);   pinMode(LEDs[2],OUTPUT);   pinMode(LEDs[3],OUTPUT);   pinMode(LEDs[4],OUTPUT);   pinMode(LEDs[5],OUTPUT);   pinMode(LEDs[6],OUTPUT);   pinMode(LEDs[7],OUTPUT);   pinMode(LEDs[8],OUTPUT);   pinMode(LEDs[9],OUTPUT);   pinMode(LEDs[10],OUTPUT);   pinMode(LEDs[11],OUTPUT);   pinMode(LEDs[12],OUTPUT);   pinMode(LEDs[13],OUTPUT);   pinMode(LEDs[14],OUTPUT);   pinMode(LEDs[15],OUTPUT);   pinMode(LEDs[16],OUTPUT);   pinMode(RPMinput,INPUT); }  void loop() {   RPMvalue=analogRead(RPMinput); //Bring RPM value in   RPMvalue=constrain(RPMvalue,0,1000); //Constrain RPM values between 0 and 1000   RPM=RPMvalue*8;  //Get the 'true' RPM value   n=RPM/500;  //Get the value of the index   remainder=RPM-n*500;  //Get the remainder of the RPM for the PWM LED   remainder=constrain(remainder,0,250); //Constrain the remainder value between 0 and 250   PWMvalue=remainder/2;  //Floored PWM value for the highest LED   Serial.print(n);      switch (n) {     case 0:       digitalWrite(LEDs[1],LOW);       digitalWrite(LEDs[2],LOW);       digitalWrite(LEDs[3],LOW);       digitalWrite(LEDs[4],LOW);       digitalWrite(LEDs[5],LOW);       digitalWrite(LEDs[6],LOW);       digitalWrite(LEDs[7],LOW);       digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       break;     case 1:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],LOW);       digitalWrite(LEDs[3],LOW);       digitalWrite(LEDs[4],LOW);       digitalWrite(LEDs[5],LOW);       digitalWrite(LEDs[6],LOW);       digitalWrite(LEDs[7],LOW);       digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       break;     case 2:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       //digitalWrite(LEDs[3],LOW);       digitalWrite(LEDs[4],LOW);       digitalWrite(LEDs[5],LOW);       digitalWrite(LEDs[6],LOW);       digitalWrite(LEDs[7],LOW);       digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[3],PWMvalue);       break;     case 3:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       //digitalWrite(LEDs[4],LOW);       digitalWrite(LEDs[5],LOW);       digitalWrite(LEDs[6],LOW);       digitalWrite(LEDs[7],LOW);       digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[4],PWMvalue);       break;     case 4:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       //digitalWrite(LEDs[5],LOW);       digitalWrite(LEDs[6],LOW);       digitalWrite(LEDs[7],LOW);       digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[5],PWMvalue);       break;     case 5:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       //digitalWrite(LEDs[6],LOW);       digitalWrite(LEDs[7],LOW);       digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[6],PWMvalue);       break;     case 6:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       //digitalWrite(LEDs[7],LOW);       digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[7],PWMvalue);       break;     case 7:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       //digitalWrite(LEDs[8],LOW);       digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[8],PWMvalue);       break;     case 8:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       //digitalWrite(LEDs[9],LOW);       digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[9],PWMvalue);       break;     case 9:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       digitalWrite(LEDs[9],HIGH);       //digitalWrite(LEDs[10],LOW);       digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[10],PWMvalue);       break;     case 10:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       digitalWrite(LEDs[9],HIGH);       digitalWrite(LEDs[10],HIGH);       //digitalWrite(LEDs[11],LOW);       digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[11],PWMvalue);       break;     case 11:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       digitalWrite(LEDs[9],HIGH);       digitalWrite(LEDs[10],HIGH);       digitalWrite(LEDs[11],HIGH);       //digitalWrite(LEDs[12],LOW);       digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[12],PWMvalue);       break;     case 12:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       digitalWrite(LEDs[9],HIGH);       digitalWrite(LEDs[10],HIGH);       digitalWrite(LEDs[11],HIGH);       digitalWrite(LEDs[12],HIGH);       //digitalWrite(LEDs[13],LOW);       digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[13],PWMvalue);       break;     case 13:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       digitalWrite(LEDs[9],HIGH);       digitalWrite(LEDs[10],HIGH);       digitalWrite(LEDs[11],HIGH);       digitalWrite(LEDs[12],HIGH);       digitalWrite(LEDs[13],HIGH);       //digitalWrite(LEDs[14],LOW);       digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[14],PWMvalue);       break;     case 14:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       digitalWrite(LEDs[9],HIGH);       digitalWrite(LEDs[10],HIGH);       digitalWrite(LEDs[11],HIGH);       digitalWrite(LEDs[12],HIGH);       digitalWrite(LEDs[13],HIGH);       digitalWrite(LEDs[14],HIGH);       //digitalWrite(LEDs[15],LOW);       digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[15],PWMvalue);       break;     case 15:       digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrite(LEDs[5],HIGH);       digitalWrite(LEDs[6],HIGH);       digitalWrite(LEDs[7],HIGH);       digitalWrite(LEDs[8],HIGH);       digitalWrite(LEDs[9],HIGH);       digitalWrite(LEDs[10],HIGH);       digitalWrite(LEDs[11],HIGH);       digitalWrite(LEDs[12],HIGH);       digitalWrite(LEDs[13],HIGH);       digitalWrite(LEDs[14],HIGH);       digitalWrite(LEDs[15],HIGH);       //digitalWrite(LEDs[16],LOW);       analogWrite(LEDs[16],PWMvalue);       break;     case 16:       for (int brightness = 0; brightness < 255; brightness++) {       analogWrite(LEDs[1],brightness);       analogWrite(LEDs[2],brightness);       analogWrite(LEDs[3],brightness);       analogWrite(LEDs[4],brightness);       analogWrite(LEDs[5],brightness);       analogWrite(LEDs[6],brightness);       analogWrite(LEDs[7],brightness);       analogWrite(LEDs[8],brightness);       analogWrite(LEDs[9],brightness);       analogWrite(LEDs[10],brightness);       analogWrite(LEDs[11],brightness);       analogWrite(LEDs[12],brightness);       analogWrite(LEDs[13],brightness);       analogWrite(LEDs[14],brightness);       analogWrite(LEDs[15],brightness);       analogWrite(LEDs[16],brightness);       delay(3);       }        for (int brightness = 255; brightness >= 0; brightness--) {       analogWrite(LEDs[1],brightness);       analogWrite(LEDs[2],brightness);       analogWrite(LEDs[3],brightness);       analogWrite(LEDs[4],brightness);       analogWrite(LEDs[5],brightness);       analogWrite(LEDs[6],brightness);       analogWrite(LEDs[7],brightness);       analogWrite(LEDs[8],brightness);       analogWrite(LEDs[9],brightness);       analogWrite(LEDs[10],brightness);       analogWrite(LEDs[11],brightness);       analogWrite(LEDs[12],brightness);       analogWrite(LEDs[13],brightness);       analogWrite(LEDs[14],brightness);       analogWrite(LEDs[15],brightness);       analogWrite(LEDs[16],brightness);       delay(3);       break;     }   } }  /* digitalWrite(LEDs[1],HIGH);       digitalWrite(LEDs[2],HIGH);       digitalWrite(LEDs[3],HIGH);       digitalWrite(LEDs[4],HIGH);       digitalWrit

Posted May 10, 2010 @ 18:14

// 0 Comments

Python saves the day!

Oh python...

Learning actual python has been challenging yet rewarding, and now its starting to pay off. I just got an Arduino and have no clue how to program it, nor do I have any idea how its math functions work. So I started pondering how to get LED's to stay on for my digital tach while allowing the highest one to PWM and change brightness. Since my Arduino Sketchpad skills were lacking (to say the least), I whipped out VIDLE and started programming in python. Five minutes later, viola--working code for the digital tach that allows me to get LED's to work the way I want when I transition to a physical representation of my idea.

# Test code for the digital tach

# By: Mike McD

# Created: 5/1/10

# Last updated: 5/2/10

# digitalredline.posterous.com

from __future__ import division # imports regular division, that is in "normal" python, 9/4 = 2, instead of 2.25. This makes it 2.25.

rpm = input("what is the desired RPM?") # input desired RPM from shell

rpmList=[0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000] # list of all possibly RPM values, in 500 RPM increments. Each corresponds to a single LED.

newList=[] # new list which stores values of LED's which should be illuminated.

n=rpm//500 # floor divison operator. Now that division works normally, I need to floor the values to get whole numbers to look things up in the list.

remainder=rpm-n*500 # this is the remainder of the division, spit out as a number between 0 and 500, the range of each LED in RPM.

PWMvalue=remainder/2 # this is the value of the Arduino PWM band, corresponding to an incriment of 0-250 in 1/2 step incriments.

newList.append(rpmList[:n+1]) # this is the list of values for LED's which should be illuminated

print rpmList[:n+1]

print newList

print remainder, PWMvalue

Posted May 2, 2010 @ 14:46

// 0 Comments

Danger, Danger!

Here is a little chunk of code I just copied and changed for my Arduino Mega board

that I bought and finally had time to start playing with.

This is literally EXACTLY what i've been wanting.

Easy PWM control makes flashing lights super easy,

so this should (in theory) make the digital tach even cooler.

Assuming I can find the right wire to tap to get the tach signal.

Anyways, this code makes an LED pulse on and off (as per those warning lights that always appear in moves).

Not on/off, but it gently pulses, warning you of impending danger.

I think I want one of these to warn me of my danger!

/* Pulsing LED, as per a warning light!  * By Mike McD  * Created on 4/30/10, Last Updated 4/30/10  */  const int ledPin = 10;  void setup() {     pinMode(ledPin, OUTPUT);}  void loop() {     // fade the LED on ledPin from off to brightest:     for (int brightness = 0; brightness < 255; brightness++){       analogWrite(ledPin, brightness);       delay(1);     }     // fade the LED on ledPin from brithstest to off:     for (int brightness = 255; brightness >= 0; brightness--) {       analogWrite(ledPin, brightness);       delay(1);     }  }

P.S. I just found out that the Arduino Sketchpad has HTML copy/paste ability.  
Finally my code gets posted in the correct format.  
I wonder if IDLE or VIDLE has this capability?
 
P.P.S Doesn't seem like it...at least not VIDLE?  Maybe some other IDE?

Posted April 30, 2010 @ 16:57

// 0 Comments

And heres the video...and a few pictures!

P.S. for those of you using macs, I used snapzprox (http://www.ambrosiasw.com/utilities/snapzprox/) to take my videos. They have like a 30 day free trial, and it does exactly what you need. Quicktime 10 (OS 10.6 and later I believe) also has a screen record feature, but since my graphics car doesn't support pixel shader 3.0 (which is necessary for the planet skins) I was unable to utilize said features.

Along with a few close ups!

Launch! And the exit burn.

Lunar Orbit (Its really NOT circular, its just weird. I'm gunna work on circular later...)

Moon "weight," showing the changing weight force acting on the rocket as it orbits!

Close up of the above transition into lunar orbit.

Here is the v vs t graph for the rocket.

Close up including launch, initial orbit (sinusoidal), burn (negative exponential), and subsequent orbit.

Posted March 31, 2010 @ 18:23

// 0 Comments

Orbit, Sweet Orbit.

from __future__ import division
from visual import *
from visual.graph import *

#Lunar Orbit Insertion Test
#By Mike McD
#Updated 3/31/10

lmlengthconstant=13

#time

dt=.5
t=0

#stage burn time
SICtime=150
SIItime=367
SIVBtime=475
SIVBtime2=335
CMburntime=3800

#total elapsed time
elSICtime=SICtime
elSIItime=SICtime+SIItime
elSIVBtime=SICtime+SIItime+SIVBtime
el1stburnstart=22000
el1stburnend=el1stburnstart+SIVBtime2
el2ndburnstart=188000
el2ndburnend=el2ndburnstart+CMburntime

#thrust
SICthrust=33400000
SIIthrust=5115000
SIVBthrust=1001000
thrust=0

#mass
SICmass=2279770.49
SIImass=480782.81
SIVBmass=119784.89
payload=47000

#drag and atmo pressure data
cd=.4
farea=78.54
scaleheight=9000
atmoinit=1.2

#earth and moon data
emdistance=384399e3

earth=sphere()
earth.radius=6.38e6
earth.opacitz=1
earth.mass=6e24
earth.material=materials.earth

moon=sphere()
moon.pos=vector(0,0,emdistance)
moon.radius=1737e3
moon.mass=7.3477e22
moon.color=(0.5,0.5,0.5)
moon.material=materials.marble

moonperiod=27.3*86400
phasedays=2.63
moonphase=phasedays*86400

moontrace=curve(color=(0.5,0.5,0.5))

#gravity constants
gconstant=6.67e-11

#initial angles and pitch program
position=0 #radians
angularposition=0 #radians

angularvelocity=0 #radians/sec

angularacceleration1=-radians(0.01) #rad/sec^2 for 250AA1starttime=130
AA1endtime=330
angularacceleration2=radians(0.007) #rad/sec^2 for 350AA2starttime=331
AA2endtime=435

finalangle=-radians(80)

#rocket frames 1-4 (Full, 2nd stage, 3rd stage, LM)
rocket=frame()
SIC=cylinder(frame=rocket,pos=(0,0,0),length=42, radius=5, color=color.white)
SII=cylinder(frame=rocket,pos=(SIC.length,0,0),length=24.9, radius=5, color=color.white)
SIV=cylinder(frame=rocket,pos=(SIC.length+SII.length,0,0),length=17.8, radius=3.3, color=color.white)
shroud=cone(frame=rocket,pos=(SIC.length+SII.length,0,0),length=10, radius=5, color=color.white)
LM=cone(frame=rocket,pos=(SIC.length+SII.length+SIV.length,0,0),length=21.5, radius=3.3,color=color.white)
SM=cylinder(frame=rocket,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant,0,0),length=7.95, radius=1.95, color=color.white)
CM=cone(frame=rocket,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant+SM.length,0,0),length=3.47, radius=1.95, color=color.blue)
engine1=cone(frame=rocket,pos=(-4,0,-3.5),length=5.79, radius=1.88, color=color.red)
engine2=cone(frame=rocket,pos=(-4,0,3.5),length=5.79, radius=1.88, color=color.red)
engine3=cone(frame=rocket,pos=(-4,0,0),length=5.79, radius=1.88, color=color.red)
engine4=cone(frame=rocket,pos=(-4,3.5,0),length=5.79, radius=1.88, color=color.red)
engine5=cone(frame=rocket,pos=(-4,-3.5,0),length=5.79, radius=1.88, color=color.red)
engine1a=cone(frame=rocket,pos=(SIC.length-2.5,0,-3.5),length=3.38, radius=1, color=color.red)
engine2a=cone(frame=rocket,pos=(SIC.length-2.5,0,3.5),length=3.38, radius=1, color=color.red)
engine3a=cone(frame=rocket,pos=(SIC.length-2.5,0,0),length=3.38, radius=1, color=color.red)
engine4a=cone(frame=rocket,pos=(SIC.length-2.5,3.5,0),length=3.38, radius=1, color=color.red)
engine5a=cone(frame=rocket,pos=(SIC.length-2.5,-3.5,0),length=3.38, radius=1, color=color.red)
engine6=cone(frame=rocket,pos=(SIC.length+SII.length-2,0,0),length=3.38, radius=1, color=color.red)
rocket.axis=(0,0,1)
rocket.pos=(0,0,earth.radius)

rocket2=frame()
SII=cylinder(frame=rocket2,pos=(SIC.length,0,0),length=24.9, radius=5, color=color.white)
SIV=cylinder(frame=rocket2,pos=(SIC.length+SII.length,0,0),length=17.8, radius=3.3, color=color.white)
shroud=cone(frame=rocket2,pos=(SIC.length+SII.length,0,0),length=10, radius=5, color=color.white)
LM=cone(frame=rocket2,pos=(SIC.length+SII.length+SIV.length,0,0),length=21.5, radius=3.3,color=color.white)
SM=cylinder(frame=rocket2,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant,0,0),length=7.95, radius=1.95, color=color.white)
CM=cone(frame=rocket2,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant+SM.length,0,0),length=3.47, radius=1.95, color=color.blue)
engine1a=cone(frame=rocket2,pos=(SIC.length-2.5,0,-2),length=3.38, radius=1, color=color.red)
engine2a=cone(frame=rocket2,pos=(SIC.length-2.5,0,2),length=3.38, radius=1, color=color.red)
engine3a=cone(frame=rocket2,pos=(SIC.length-2.5,0,0),length=3.38, radius=1, color=color.red)
engine4a=cone(frame=rocket2,pos=(SIC.length-2.5,2,0),length=3.38, radius=1, color=color.red)
engine5a=cone(frame=rocket2,pos=(SIC.length-2.5,-2,0),length=3.38, radius=1, color=color.red)
engine6=cone(frame=rocket2,pos=(SIC.length+SII.length-2,0,0),length=3.38, radius=1, color=color.red)
rocket2.axis=(0,0,1)
rocket2.pos=(0,0,earth.radius)

rocket3=frame()
SIV=cylinder(frame=rocket3,pos=(SIC.length+SII.length,0,0),length=17.8, radius=3.3, color=color.white)
shroud=cone(frame=rocket3,pos=(SIC.length+SII.length,0,0),length=10, radius=5, color=color.white)
LM=cone(frame=rocket3,pos=(SIC.length+SII.length+SIV.length,0,0),length=21.5, radius=3.3,color=color.white)
SM=cylinder(frame=rocket3,pos=(SIC.length+SII.length+SIV.length-lmlengthconstant+LM.length-lmlengthconstant,0,0),length=7.95, radius=1.95, color=color.white)
CM=cone(frame=rocket3,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant+SM.length,0,0),length=3.47, radius=1.95, color=color.blue)
engine6=cone(frame=rocket3,pos=(SIC.length+SII.length-2,0,0),length=3.38, radius=1, color=color.red)
rocket3.axis=(0,0,1)
rocket3.pos=(0,0,earth.radius)

rocket4=frame()
LM=cone(frame=rocket4,pos=(SIC.length+SII.length+SIV.length,0,0),length=21.5, radius=3.3,color=color.white)
SM=cylinder(frame=rocket4,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant,0,0),length=7.95, radius=1.95, color=color.white)
CM=cone(frame=rocket4,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant+SM.length,0,0),length=3.47, radius=1.95, color=color.blue)
engine8=cone(frame=rocket4,pos=(SIC.length+SII.length+SIV.length+LM.length-lmlengthconstant-5,0,0),length=3.82, radius=1.25, color=color.red)
rocket4.axis=(0,0,1)
rocket4.pos=(0,0,earth.radius)

#trace and initial positions
trace=curve(color=color.orange)
trace2=curve(color=color.cyan)
trace3=curve(color=color.magenta)
trace4=curve(color=color.white)
trace5=curve(color=color.blue)
trace6=curve(color=color.red)
rocket.pos=(0,0,earth.radius)
rocket2.pos=(0,0,earth.radius)
rocket3.pos=(0,0,earth.radius)
rocket4.pos=(0,0,earth.radius)
rocket.velocity=vector(464,0,0)

#wind velocity
boost=(2*pi*(sqrt((rocket.pos.x)**2)+((rocket.pos.z)**2)))/(86400)

#altitude
altitude1=sqrt(((rocket.pos.x)**2)+((rocket.pos.z)**2))-earth.radius
altitude2=sqrt(((rocket2.pos.x)**2)+((rocket2.pos.z)**2))-earth.radius
altitude3=sqrt(((rocket3.pos.x)**2)+((rocket3.pos.z)**2))-earth.radius
altitude4=sqrt(((rocket4.pos.x)**2)+((rocket4.pos.z)**2))-earth.radius

#scene data
scene=display(center=(rocket4.pos.x,rocket4.pos.z,rocket4.pos.z))
scene.select()

#display 1
graph1=gdisplay(x=0,y=0,width=600,height=400,title='moongrav vs time',xtitle='time',ytitle='moongrav',xmax=225000,xmin=0,ymax=7.5e3,ymin=0,foreground=color.white,background=color.black)
rocketposangle=gcurve(gdisplay=graph1,color=color.white)

#display2
graph1=gdisplay(x=0,y=0,width=600,height=400,title='velocity vs time',xtitle='time',ytitle='velocity',xmax=200000,xmin=190000,ymax=20000,ymin=0,foreground=color.white,background=color.black)
velocity=gcurve(gdisplay=graph1,color=color.white)

while true:
rate=2000
#moon movement
t=t+dt

#if tAA2endtime:
angularacceleration=0
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=finalangle
else:
angularacceleration=0
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=finalangle
#rocket movement
if t #time
t=t+dt
#angles
position=atan2(rocket.pos.z,rocket.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket.pos)
rocket2.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket2.pos)
rocket3.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket3.pos)
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket.pos.x)**2)+((rocket.pos.z)**2)))/(86400)
boostx=boost*sin(position)
boostz=boost*cos(position)
#mass and weight
SICmass=SICmass*(1-(t/SICtime))
SIImass=SIImass
SIVBmass=SIVBmass
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=((gconstant*masstotal*earth.mass)/(((rocket.pos.x)**2)+((rocket.pos.z)**2)))/(masstotal)
weight=masstotal*grav*-1
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=SICthrust
thrustx=thrust*cos(absolute)
thrustz=thrust*sin(absolute)
#drag
air=atmoinit*(e**((-(altitude1)/(scaleheight))))
drag=.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon)**2)+((ztomoon)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#movement
trace.append(pos=rocket.pos, retain=100000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude1=sqrt(((rocket.pos.x)**2)+((rocket.pos.z)**2))-earth.radius
rocket.pos=rocket.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket2.pos=rocket2.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket3.pos=rocket3.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon,xtomoon)
#print position
#print t
#print degrees(absolute)
#print boostx
#print boostz
#print rtomoon
else:
if elSICtime #time
t=t+dt
#angles
position=atan2(rocket2.pos.z,rocket2.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket2.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket2.pos)
rocket3.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket3.pos)
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket2.pos.x)**2)+((rocket2.pos.z)**2)))/(86400)
boostx=boost*sin(position)
boostz=boost*cos(position)
#mass and weight
SICmass=0
SIImass=SIImass*(1-((t-SICtime)/(SIItime)))
SIVBmass=SIVBmass
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=((gconstant*masstotal*earth.mass)/(((rocket2.pos.x)**2)+((rocket2.pos.z)**2)))/(masstotal)
weight=masstotal*grav*-1
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=SIIthrust
thrustx=thrust*cos(absolute)
thrustz=thrust*sin(absolute)
#drag
air=atmoinit*(e**((-(altitude2)/(scaleheight))))
drag=.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon2)**2)+((ztomoon2)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#movement
trace2.append(pos=rocket2.pos, retain=10000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude2=sqrt(((rocket2.pos.x)**2)+((rocket2.pos.z)**2))-earth.radius
rocket2.pos=rocket2.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket3.pos=rocket3.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket.visible=false
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon2,xtomoon2)
#print position
#print t
#print degrees(absolute)
#print boostx
#print boostz
#print rtomoon
else:
if elSIItime #time
t=t+dt
#angles
position=atan2(rocket3.pos.z,rocket3.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket3.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket3.pos)
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket3.pos.x)**2)+((rocket3.pos.z)**2)))/(86400)
boostx=boost*sin(position)
boostz=boost*cos(position)
#mass and weight
SICmass=0
SIImass=0
SIVBmass=SIVBmass*(1-((t-SICtime-SIItime)/(SIVBtime)))
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=-1*((gconstant*masstotal*earth.mass)/(((rocket3.pos.x)**2)+((rocket3.pos.z)**2)))/(masstotal)
weight=masstotal*grav
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=SIVBthrust
thrustx=thrust*cos(absolute)
thrustz=thrust*sin(absolute)
#drag
air=atmoinit*(e**((-(altitude3)/(scaleheight))))
drag=.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon3)**2)+((ztomoon3)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#movement
trace3.append(pos=rocket3.pos, retain=10000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude3=sqrt(((rocket3.pos.x)**2)+((rocket3.pos.z)**2))-earth.radius
rocket3.pos=rocket3.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket2.visible=false
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon3,xtomoon3)
#print position
#print t
#print degrees(absolute)
#print boostx
#print boostz
else:
if elSIVBtime #time
t=t+dt
#angles
position=atan2(rocket4.pos.z,rocket4.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(86400)
boostx=boost*cos(position)
boostz=boost*sin(position)
#mass and weight
SICmass=0
SIImass=0
SIVBmass=0
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=-1*((gconstant*masstotal*earth.mass)/(((rocket.pos.x)**2)+((rocket.pos.z)**2)))/(masstotal)
weight=masstotal*grav
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=0
thrustx=thrust*cos(absolute)
thrustz=thrust*sin(absolute)
#drag
#air=atmoinit*(e**((-(altitude4)/(scaleheight))))
drag=0 #.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon4)**2)+((ztomoon4)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#motion
trace4.append(pos=rocket4.pos, retain=25000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude4=sqrt(((rocket4.pos.x)**2)+((rocket4.pos.z)**2))-earth.radius
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket3.visible=false
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon4,xtomoon4)
#print position
#print t
#print degrees(absolute)
#print boostx
#print boostz
else:
if el1stburnstart #time
t=t+dt
#angles
position=atan2(rocket4.pos.z,rocket4.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(86400)
boostx=boost*cos(position)
boostz=boost*sin(position)
#mass and weight
SICmass=0
SIImass=0
SIVBmass=0
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=-1*((gconstant*masstotal*earth.mass)/(((rocket.pos.x)**2)+((rocket.pos.z)**2)))/(masstotal)
weight=masstotal*grav
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=SIVBthrust
thrustx=thrust*cos(absolute)
thrustz=thrust*sin(absolute)
#drag
#air=atmoinit*(e**((-(altitude4)/(scaleheight))))
drag=0 #.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon4)**2)+((ztomoon4)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#motion
trace5.append(pos=rocket4.pos, retain=25000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude4=sqrt(((rocket4.pos.x)**2)+((rocket4.pos.z)**2))-earth.radius
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket3.visible=false
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon4,xtomoon4)
#print position
#print t
#print degrees(absolute)
#print boostx
#print boostz
else:
if el1stburnend #time
t=t+dt
#angles
position=atan2(rocket4.pos.z,rocket4.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(86400)
boostx=boost*cos(position)
boostz=boost*sin(position)
#mass and weight
SICmass=0
SIImass=0
SIVBmass=0
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=-1*((gconstant*masstotal*earth.mass)/(((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(masstotal)
weight=masstotal*grav
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=0
thrustx=thrust*cos(absolute)
thrustz=thrust*sin(absolute)
#drag
#air=atmoinit*(e**((-(altitude4)/(scaleheight))))
drag=0 #.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon4)**2)+((ztomoon4)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#motion
trace4.append(pos=rocket4.pos, retain=25000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude4=sqrt(((rocket4.pos.x)**2)+((rocket4.pos.z)**2))-earth.radius
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket3.visible=false
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon4,xtomoon4)
#print position
#print t
#print degrees(absolute)
#print boostx
#print boostz
else:
if el2ndburnstart #time
t=t+dt
#angles
position=atan2(rocket4.pos.z,rocket4.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(86400)
boostx=boost*cos(position)
boostz=boost*sin(position)
#mass and weight
SICmass=0
SIImass=0
SIVBmass=0
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=-1*((gconstant*masstotal*earth.mass)/(((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(masstotal)
weight=masstotal*grav
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=100000
thrustx=thrust*cos(absolute-1.8)
thrustz=thrust*sin(absolute-1.8)
#drag
#air=atmoinit*(e**((-(altitude4)/(scaleheight))))
drag=0 #.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon4)**2)+((ztomoon4)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#motion
trace6.append(pos=rocket4.pos, retain=25000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude4=sqrt(((rocket4.pos.x)**2)+((rocket4.pos.z)**2))-earth.radius
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket3.visible=false
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon4,xtomoon4)
#print degrees(position)
#print t
#print moonweightx
#print moonweightz
#print degrees(absolute)
#print rtomoon4
#print boostx
#print boostz
else:
#time
t=t+dt
#angles
position=atan2(rocket4.pos.z,rocket4.pos.x)
angularvelocity=angularvelocity+angularacceleration*dt
angularposition=angularposition+angularvelocity*dt+angularacceleration*dt*dt
absolute=position+angularposition
#rotation
rocket4.rotate(angle=angularvelocity*dt,axis=(0,1,0),origin=rocket4.pos)
#wind velocity
boost=0 #(2*pi*(sqrt((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(86400)
boostx=boost*cos(position)
boostz=boost*sin(position)
#mass and weight
SICmass=0
SIImass=0
SIVBmass=0
payload=payload
masstotal=SICmass+SIImass+SIVBmass+payload
grav=-1*((gconstant*masstotal*earth.mass)/(((rocket4.pos.x)**2)+((rocket4.pos.z)**2)))/(masstotal)
weight=masstotal*grav
weightx=weight*cos(position)
weightz=weight*sin(position)
#thrust
thrust=0
thrustx=thrust*cos(absolute-1.8)
thrustz=thrust*sin(absolute-1.8)
#drag
#air=atmoinit*(e**((-(altitude4)/(scaleheight))))
drag=0 #.5*air*(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))*farea*cd
dragx=drag*cos(absolute)
dragz=drag*sin(absolute)
#moon gravity
moongrav=((gconstant*masstotal*moon.mass)/(((xtomoon4)**2)+((ztomoon4)**2)))/(masstotal)
moonweightx=masstotal*moongrav*cos(moontheta)
moonweightz=masstotal*moongrav*sin(moontheta)
#net force x
netforcex=thrustx+weightx-dragx+moonweightx
#net force z
netforcez=thrustz+weightz-dragz+moonweightz
#net acceleration x
netaccelerationx=netforcex/masstotal
#net accleration z
netaccelerationz=netforcez/masstotal
#motion
trace4.append(pos=rocket4.pos, retain=2500000)
rocket.velocity.x=rocket.velocity.x+netaccelerationx*dt
rocket.velocity.z=rocket.velocity.z+netaccelerationz*dt
altitude4=sqrt(((rocket4.pos.x)**2)+((rocket4.pos.z)**2))-earth.radius
rocket4.pos=rocket4.pos+(rocket.velocity+(boostx,0,boostz))*dt
rocket3.visible=false
#plot
rocketposangle.plot(pos=(t,(sqrt(((moonweightx)**2)+((moonweightz)**2)))))
velocity.plot(pos=(t,(sqrt(((rocket.velocity.x)**2)+((rocket.velocity.z)**2)))))
moontheta=atan2(ztomoon4,xtomoon4)
#print degrees(position)
#print t
#print moonweightx
#print moonweightz
#print degrees(absolute)
#print rtomoon4
#print boostx
#print boostz

SWEEEEEEET! Finally! This isn't exactly the stable circular orbit that apollo astronauts wanted, but it does work. I will work on a circular orbit on the multitude of plane flights that I have over break.

Overview of the project:
Re-create the lunar lander project from excel into visual python, due to heavy prompting from peers (It's just a graph? Nothing moves? LAME)
Build a scale model of the rocket in VP, use the force prospective rather than momentum to make a more realistic view of the launch and subsequent orbit.
Get it into a stable lunar orbit!

I will post some more pictures later, along with a video!

Posted March 31, 2010 @ 12:53

// 0 Comments

Physical Pendulum!

Here is the code for my physical pendulum.

from __future__ import division # makes sure is decimal and not integer
from visual import *
from visual.graph import * # graphing capability

Cross=frame()
Yardstick=box(frame=Cross,pos=(0,0,0),material=materials.wood,size=(.028,1,.008))
Crosspiece=cylinder(frame=Cross,pos=(-.267,-.3,.009525),material=materials.plastic,radius=.009525,length=.534)
BattL=box(frame=Cross,pos=(-.275,-.3,.009525),size=(.016,.043,.025),color=color.blue)
BattR=box(frame=Cross,pos=(.275,-.3,.009525),size=(.016,.043,.025),color=color.blue)
Cross.origin=(0,0,0)

#starting positionstartAngle = 20
theta = math.radians(startAngle)
Cross.rotate(angle = theta, axis = (0,0,1), origin = (0,.5,0))
angVel = 0
angAcc = 0

# masses
Yardstick.mass=0.1637755102
Crosspiece.mass=0.133877551
BattL.mass=0.0478061224
BattR.mass=0.0478061224
massTotal=Yardstick.mass+Crosspiece.mass+BattL.mass+BattR.mass

# moments
cmLength=(Yardstick.mass*(Yardstick.size.y/2)+((Crosspiece.mass+BattL.mass+BattR.mass)*(0.8*Yardstick.size.y)))/(Yardstick.mass+Crosspiece.mass+BattL.mass+BattR.mass)
rotInertia=((1/3)*Yardstick.mass*Yardstick.length)+((1/12)*Crosspiece.mass*((0.8)**2))+(2*((0.275)**2))+((Crosspiece.mass+BattL.mass+BattR.mass)*((0.8)**2))

# set up angle graph
graph1 = gdisplay(x=150, y=600, width=400, height=300,
title='Angle vs. Time', xtitle='time (s)', ytitle='angle (rad)',
xmax=10., xmin=0., ymax=1.1*theta, ymin=-1.1*theta,
foreground=color.black, background=color.white)
acurve = gcurve(gdisplay = graph1, color = color.red)

# time
time = 0
dt = 0.01

while time <4*(2*pi*sqrt(rotInertia/(massTotal*9.8*cmLength))):
rate(100)
time += dt
gravTorque = -1 * cmLength * (Yardstick.mass+Crosspiece.mass+BattL.mass+BattR.mass) * 9.8 * sin(theta)
angAcc = gravTorque / rotInertia
angVel += angAcc * dt
theta += angVel * dt
Cross.rotate(angle = angVel * dt, axis = (0,0,1), origin = (0,.5,0))
acurve.plot(pos = (time, theta))
#print rotInertia
#print 2*pi*sqrt(rotInertia/(massTotal*9.8*cmLength))

Here is a picture of the "pendulum" cross thing:

Here is the cosine graph of the function!

Fun facts about the pendulum: The center of mass is located 0.675 meters from the tip of the yard stick, or just above where the cross is, and a moment of inertia of 0.307 kg*m^2. Both of these factors lead to a period of 2.158 seconds.

Posted March 25, 2010 @ 16:16

// 0 Comments

Shoot the moon.

from __future__ import division
from visual import *
from visual.graph import *

#Vertical Launch Test
#By Mike McD
#Updated 3/22/10

lmlengthconstant=13

#time
dt=.5
t=0

#stage burn time
SICtime=150
SIItime=367
SIVBtime=475
SIVBtime2=335
CMburntime=500

#total elapsed time
elSICtime=SICtime
elSIItime=SICtime+SIItime
elSIVBtime=SICtime+SIItime+SIVBtime
el1stburnstart=22000
el1stburnend=el1stburnstart+SIVBtime2
el2ndburnstart=190200
el2ndburnend=el2ndburnstart+CMburntime

#thrust
SICthrust=33400000
SIIthrust=5115000
SIVBthrust=1001000
thrust=0

#mass
SICmass=2279770.49
SIImass=480782.81
SIVBmass=119784.89
pazload=47000

#drag and atmo pressure data
cd=.4
farea=78.54
scaleheight=9000
atmoinit=1.2

#earth and moon data
emdistance=384399e3

earth=sphere()
earth.radius=6.38e6
#earth.color=color.blue
earth.opacitz=1
earth.mass=6e24
earth.material=materials.earth

moon=sphere()
moon.pos=vector(emdistance,0,0)
moon.radius=1737e3
moon.mass=7.3477e22
moon.color=color.yellow

moonperiod=27.3*86400
phasedays=2.55
moonphase=phasedays*86400

moontrace=curve(color=color.yellow)

#gravity constants
gconstant=6.67e-11

#initial angles and pitch program
position=0 #radians
angularposition=0 #radians

angularvelocity=0 #radians/sec

angularacceleration1=-radians(0.01) #rad/sec^2 for 250AA1starttime=130
AA1endtime=330
angularacceleration2=radians(0.007) #rad/sec^2 for 350AA2starttime=331
AA2endtime=435

finalangle=-radians(80)
TLIangle=-radians(36.12+180)