/* Possible Changes :
   -> temp variables in functions?
   -> redundant distrib/distribution syntax !
   -> UI ability to change label - element, component, unit, etc
   -> Function overloading?
   -> When it is 'key mass', the origin should default to the current component,
      not root
   -> add keyword parent, as a reference, eg key mass of parent
*/

default mass = real 0;
default power = real 0;
default cost = real 0;
default efficiency = real 1;

/* #include "general.c" */
Pi = 3.14159265358979323846;
G  = 9.81; /* m/s^2 */

Component LunarIce {
  
  Component Mission /* Baseline */
    {
      sortieDuration = 20; /* hours */
      duration = 30*24;    /* hours */
      distrib speed = distribution(min=0.20, max=0.4, avg=0.25);
      distance = duration*3600*get_min(Speed)*2;
      Component Launcher
	{
	  FairingDia = 2.3; /* 2.3 meter payload fairing */
	  InjectionError = 1.5; /* deg */
	}
    }
    
  component Earth 
    {
      /* this would normally be added by the robot design rather
	 than in the component itself. Do when we have ability
	 to reference & include components */
      type CommLink; 

      g = 9.81;
      day = 24; /* hours */

      Component Atmosphere { CommAmplitude = -0.3; }
      Component Ionosphere { CommAmplitude = -0.1; }
      Component Rain {
	distrib CommAmplitude = distribution(max=0, min=-1, avg=-0.2); }
    }

  component Moon 
    {
      RevPeriod = 28.45 * 24 * 3600; /* in seconds */
      distrib EarthDist = Distribution(avg=384.0E6, max=406.697E6, min=356.410E6); 
      /* meters */
      
      day = 27.322 + 10/24; /* approx max "day" hrs - sun cycle, see p57 LSB */
      g = 9.81/6;
    }
  
  component Space 
    {
      Temp = 3; /* degree Kelvin */
    }

  Component Spacecraft {
    transparent mass = sum(key mass of Spacecraft) : 1500;
    
    Component Icebreaker {
	
      transparent mass = sum(key mass of icebreaker) : 392; /* kg */
      transparent power = sum(key power of icebreaker) : 400; /* W */

      Component Locomotion {
	Component Spec 
	  {
	    slopemax = 22 degree; 
	    obstacleMax = 0.25;   /* meters  */
	    turnmin = 6.0; /* min outer turning radius */
	  }

	Component Soil 
	  { /* Lunar Soil LSS model, LSB, P526,8,.. */
	    /* Pressure sinkage characterstics under wheel load */
	    n = 1; /* regolith sinkage exponent */
	    Kc = 1382.32; /* cohesive modulus N/m^(n+1) */
	    Kphi = 816331.27; /* frictional modulus N/m^(n+2) */
	    
	    c = 172.79; /* cohesion of regolith N/m^2 */
	    K = 0.01778; /* coeff of slip  meters */
	    phi = 35*Pi/180; /* angle of internal friction/ repose 35 */
	    fr  = 0.1; /* coeffecient of rolling friction */
	    pressureMax = 6000; /* Pascals */
	  }

	CGz = 1.0; /* meters offset from bottom center */
	CGL = 0.1; /* length & width wise horizontal offsets */
	CGw = 0.1;
	
	L = fairingDia/sqrt(2) - 0.01; /* length meters */
	W = sqrt(fairingDia^2 - (L+0.01)^2) - 0.01; /* width meters */
	n = 4; /* n wheels */
	
	Contact(load) = sqrt(wheel.d^2/4.0 - 
			     (wheel.d/2 - sinkage(load)^2));
	contactArea(load) = Contact(load)*wheel.w;
	pressureMax = (old icebreaker.mass)*moon.G/
	  (2*contactArea(WheelLoadMax));
	slopeLmax = atan((L/2 - wheel.d/2 - CGL)/CGz);
	slopeWmax = atan((W/2 - CGw)/CGz);
	slopemax = max(slopeLmax, slopeWmax);

	weight = (old icebreaker.mass)*moon.G;
	load(theta=0,w=weight) = w/n*cos(theta)*1.1; 
	/* 1.1 accounts for imbalance, cross slopes */
	wheelLoadMax = weight/2;
	  

	/* Total electric Driving power in W */

	P(theta=0,V=mission.speed) = amplifier.P(theta,V)*n;

	/* Now calc electric power in for different conditions ??? */
	/* Drive sizing to include flat and slope conditions ??? */

	mass = 0.30/(1-0.30)*(old icebreaker.mass);
	power = P(theta=0,V=get_avg(mission.speed))*1.05; /* 1.05 is margin */
	cost = 6;

	/* May be worth checking out flexible wheel for power/ size */
	Component Wheel{ /* rigid wheel model */
	  d = 0.6; /* m */
	  w = 0.3; /* m */
	  
	  /* for rigid wheel, sinkage equation */
	  sinkage(t,w=weight) = 
	    (3*load(t,w)/((3-soil.n)*(soil.kc+W^2*soil.kphi)
				   *sqrt(D)))^(2/(2*soil.n+1));
	  
	  /* Soil thrust */
	  H(t,w=weight) = (soil.c*contactArea(load(t,w)) +
			   load(t,w)*tan(soil.phi));
	  
	  /* Compaction resistance */
	  Rc(t,w) = (3*load(t,w)/sqrt(D))^
	    ((2*soil.n+2)/(2*soil.n+1))/
	      ((3-soil.n)^((2*soil.n+2)/(2*soil.n+1))*(soil.n+1)
	       *(soil.kc+W*soil.kphi)^(1/(2*soil.n+1)));
	  
	  /* Buldozing resistance */
	  alpha = 1.0;
	  Rb(t,w) = 0.5*alpha*W*sinkage(0,load(t,w))^2*tan(pi/4 + 
							     soil.phi/2)^2
	    + 2*soil.c*W*sinkage(0,load(t,w))*tan(pi/4+soil.phi/2);
	  
	  /* Rolling resistance */
	  Rr(t,w) = soil.fr*load(t,w);
	  
	  /* Gravitational resistance */
	  Rg(t,w=weight) = w/n * sin(t);
	  
	  /* Acceleration resistance */
	  Ra(a=0) = (old icebreaker.mass)/n * a;
	  
	  /* Obstacle climbing resistance */
	  Ro(size=0) = 0;
	  
	  /* Total resistance excluding obstacle, steer, acceleration */
	  Rall(t=0,w=weight) = Rc(t,w)+Rb(t,w)+Rr(t,w)+Rg(t,w);
	  
	  /* Drawbar Pull */
	  DP(t=0,w=weight) = H(t,w) - Rall(t,w);
	  
	  /* Driving Torque */
	  
	  defl = 0.0; /* deflection of wheel */
	  Tw(slope=0,w=weight) = Rall(slope,w)*(D/2-defl);
	  
	  Tmax = Tw(0,wheelLoadMax/load(w=1));
	  
	  NCycles = mission.distance/(Pi*D);

	  /* Wheel Design Constraints */
	  sinkage(WheelLoadMax) <= 0.1*wheel.D;
	  constraint sinking = pressureMax <= soil.pressureMax;
	} /* Wheel */
			     
	Component Axle
	  {type wheel;
	   efficiency = 0.8;}
	Component Bearing{type wheel;
			  efficiency = 0.9;}
	Component SteerLink{}
	
						 
	Component Amplifier
	  {
	    /* electrical power input to amplifier */
	    P(theta,V) = Wheel.Tw(theta)*
	      wheeldrive.omega(V)/efficiency/wheeldrive.Nu;
	    
	    /* At peak torque condition */
	    PTmax = wheel.Tmax*wheeldrive.omegapeak/efficiency/Wheeldrive.Nu;
	    
	    efficiency = 0.85;
	  }
	Component WheelDrive
	  {
	    inherited type wheel;
	    Component reduction1 
	      {
		efficiency = 0.7;
		ReductionRatio = 100; /* Reduction ratio */
	      }
	    Component reduction2
	      {
		efficiency = 0.7;
		ReductionRatio = 100; /* Reduction ratio */
	      }
	    safetyFactor = 1.5; /* for specing torque */
	    
	    /* Overall effeciency of wheel & drive train in delivering
	       power from motor output to terrain */
	    Nu = product(key efficiency of type wheel of locomotion)/
	      safetyFactor;

	    /* Nominal drive train spec at output */
	    Torque = wheel.Tw; 
	    omega(V=mission.speed) = V/(wheel.D/2-wheel.defl);

	    /* Peak drive train spec at output */
	    TorquePeak = wheel.Tmax;
	    omegaPeak = get_min(mission.speed)/(wheel.D/2-wheel.defl);
	    
	    MotorNcycles = 10E6; /* real numbers needed */
	    MotorNCycles > wheel.Ncycles*2.0*ReductionRatio;
	  }

	Component SteerDrive
	  {
	    efficiency=0.3;
	  }

	sqrt(L^2 + W^2) <= fairingDia;
	2.2 * wheel.D <= L;
	slopemax >= spec.slopemax;

	/* the step/D ratio is really a function of (L+-CGL)/L, 
	   traction/weight & L/D, get something better for the hack */
	D >= obstacleMax/0.25; /* What is this factor really? */
      } /* locomotion */

      Component Structure{
	mass = 0.2*(old icebreaker.mass
		       /* except desc of locomotion ??? BUG??? */ 
	  - locomotion.mass) /* HACK ??? */ ;
	cost = 2;
	Component Body{}
	Component Mast{ mass = 0.6*sum(key mass of type mastpayload);}
	Component LanderSRMInterface {
	  mass = 0.005*((old icebreaker.mass)+lander.mass);
	  cost = 0.2;
	}
	Component LanderRoverinterface {mass = 0.005*(old icebreaker.mass);}
	Component SRMLauncherInterface {mass = 0.005*(old spacecraft.mass);}
      } /* Structure */
      
      Component PowerSystem{
	Component SolarArray { /* 4 mil cell, 3 mil coating */
	  W = 2.5; 
	  H = 1.25;
	  area = H * W; /* w x h */
	  mass = 0.75 * area * 2; /* 0.75 kg/sq.m, 2 - structural overhead */
	  unitcost = 800; /* 800 dollars/watt */
	  efficiency = 0.14; 
	  powerOut = area*1358*efficiency;
	  cost = powerOut*unitcost + 0.5; /* .5 for DD&T */
	}
	Component Actuator{ /* Solar array 1 axis pointer */
	  mass = 5;
	  distrib power = distribution(max=15,avg=5,min=0);
	  cost = 0.5;
	}
	Component BatteryPack{ /* HR190 - Yardney */
	  unitMass = 1/112; /* kg/Whr */
	  unitVolume = 0.256834; /* cc/Whr */
	
	  energy = sortieDuration * old icebreaker.power;
	  mass = energy*unitMass;
	  volume = energy*unitVolume; /* cc */
	  cost = 2; /* includes engineering, DD&T, proto */
	}
	Component PMAD{ /* Power Management And Distribution */
	  cost = 1;
	  margin = 0.2; /* power margin for the overall system */
	  Component Charger{mass = 4;}
	  Component Switch{mass = 1;}
	  Component DCDC{
	    efficiency = 0.85; 
	    mass = 2;
	    dcpowerneed = old icebreaker.power - old locomotion.power
	      - old power;
	    /* actually, need to seperate dc-dc power from margin, num
	       correct though ?? */
	    power = dcpowerneed/efficiency * (1+margin) -  
	      dcpowerneed + locomotion.power*margin : 50;
	  }
	  Component Harness{mass = 3;}
	}/* PMAD */
      } /* PowerSystem */

      Component Thermal{
	cost = 2;
	Component EBox{}
	MLImass(area,layers) = area*layers*0.04; /* 1 mil (25 um - 1600 kg/m^3) kapton */
	Component MLIblanket{area = 5; layers = 10; mass = MLImass(area,layers);}
	Component MLIplumbing{area = 1; layers = 10; mass = MLImass(area,layers);}
	Component HeatPipe{mass = 2;}
	Component Radiator{mass = 3;}
	Component Heaters{mass = 0.5;}
      } /* Thermal */
      
      Component Computronics{
	Component Cage{mass = 4; cost = 0.3;}
	
	/* once there ability to spawn or prototype boards,
	   spawn N boards, and tweak power on each to on/off ??? */
	Component Boards{
	  
	  mipsGTotal = 20; /* General purpose needed */
	  mipsOTotal = 140; /* Other mips, DSP or otherwise */
	  efficiency = 0.7;
	  int boardsActive = min(mipsGTotal/bmipsG/efficiencyG+1, 
			  mipsOTotal/bmipsDSP/efficiencyDSP+1);

	  int N = boardsActive+2;
	  mass = bmass*N; /* including heat sink, 5 of these */
	  power = bpower*boardsActive;
	  cost = 1.5 + N*0.15;
	  Component Board {
	    bmass = 1; /* including heat sink */
	    bpower = 15;
	    bmipsG = 20;
	    efficiencyG = 0.5;
	    efficiencyDSP = 0.9;
	    bmipsDSP = 25*2;
	  }
	}
	Component ConnectorBus{
	  cost = 0.1;
	  mass = 0.02*200; /* about 4kg for 200 */
	}
	Component DataRecorder{
	  cost = 0.8;
	  mass = 2;
	  distrib power = distribution(max=15,avg=4,min=1);
	}
	Component Filters{ /* Analog filters & signal conditioners */
	  cost = 0.1;
	  mass = 1;
	}
	Component SafetyCircuit{
	  cost = 0.2;
	  mass = 1;
	}
      } /* Computronics */

      Component Comm
	{ /* If there is no satellite */
	  Component transciever
	    {
	      mass = 5;
	      distrib power = distribution(max=40,avg=40,min=0);
	    }
	  Component Antenna
	    { 
	      type mastpayload;
	      mass = 1;
	    }
	  Component Cable
	    {
	      length = 5; /* meters */
	      mass = length*0.5;
	    }
	  Component backup
	    {
	      mass = 3;
	      distrib power = distribution(max=30,avg=0,min=0);
	    }
	} /* Comm(Earth) */

      Component Sensors{
	Component StarTracker{
	  mass = 0.6*2;
	  power = 3.25*2;
	  cost = 1.25;
	}
	Component IMU{
	  mass = 0.7;
	  power = 5;
	  cost = 0.25;
	}
	Component CameraPanosphere{
	  mass = 3;
	  power = 15;
	  cost = 2;
	}
	Component SunTracker {}
	Component EarthTracker {}
	Component TerrainRadar {
	  driveFOV = 120; /* what do i need to drive */
	  int Nfront = driveFOV/fov*2;
	  /* *2 for double coverage */
	  Nside = (180-driveFOV)/30; /* for single coverage on sides */
	  N = Nfront*2 + Nside*2; /* *2 for front & back */
	  
	  H = 0.8; /* meters */
	  Lookahead1 = H; /* Min lookahead */
	  LookAngle1 = atan(Lookahead1/H);
	  LookAngle2 = LookAngle1 + FOV;
	  Lookahead2 = tan(lookAngle2)*H;
	  
	  Xres = get_avg(mission.speed)/freq*4; /* X is along motion */
	  Yres = Xres;
	  
	  mass  = 0.4*N;
	  power = 0.8*(Nfront + Nside*2);
	  cost  = 2;

	  Component RadarUnit
	    {
	      freq = 20; /* Hz */
	      FOV = 30; /* deg */
	    }
	}
	Component internalSensors{
	  mass = 2;
	  power = 4;
	  cost = 0.5;
	}
      } /* sensors */

      Component Software {
	Component OS {}
	Component Safeguarding 
	  {
	  }
	Component Messaging {}
      }
      
      Component SciencePackage {
	Component CryoDrill {
	  cost = 2;
	  mass = 4;
	  distrib power = distribution(max=15,avg=0,min=0);
	  Component IRFiber {}
	  Component IRTransciever {}
	}
	Component Microscope {}
	Component REGA {
	  cost = 0.5;
	  mass = 4;
	  distrib power = distribution(max=20,avg=0,min=0);
	}
	Component PanTilt {
	  cost = 1;
	  mass = 3;
	  distrib power = distribution(max=10,avg=2,min=0);
	}
	Component IRCamera {
	  cost = 0.2;
	  mass = 2;
	  distrib power = distribution(max=5,avg=2,min=0);
	}
	Component IRLaser {
	  cost = 0.2;
	  mass = 2;
	  distrib power = distribution(max=10,avg=2,min=0);
	}
	Component HiResCamera {
	  cost = 0.5;
	  mass = 1.5;
	  distrib power = distribution(max=8,avg=5,min=0);
	}
	Component Strobe {
	  cost = 0.5;
	  mass = 1;
	  distrib power = distribution(max=10,avg=4,min=0);
	}
      } /* SciencePackage */
      
    } /* IceBreaker */

  Component Lander{
    transparent mass = sum( key mass of lander ) : 200;
    cost = 15;
    Component Fuel
      {
	Isp = 220;
	dVmidcourse = 50; /* m/s, for a 1.5 deg max error, confirm */
	Mmidcourse = old spacecraft.mass - old dMmidcourse;
	dMmidcourse = rocket.DeltaM(Mmidcourse, Isp, dVmidcourse) : 35;
	
	// dVacs = ; /* m/s, attitude control during SRM burn */
	Macs = Mmidcourse - old dMacs;
	dMacs = Mmidcourse/1456*25; /* kg */

	dVterminal = 200; /* m/s */
	Mterminal = Macs - SRM.mass - old dMterminal;
	dMterminal = rocket.DeltaM(Mterminal, Isp, dVterminal);

	Mhover = Mterminal + dMterminal;
	hangtime = 90; /* seconds of course */
	mdotfrac = Isp*Earth.G/Moon.G; 
	/* ratio of fuel mass/s to craft mass */
	dMhover = hangtime*mdotfrac/Mterminal : 35; 
	/* kg, for hovering above surface */

	dMmargin = 0.20*mass : 25; /* kg, fuel margin */

	mass = sum (key "dM*" of fuel) : 130; /* kg hydrazine fuel */
      }
    Component Propulsion
      {
	Component DE 
	  {
	    mass = 60; /* Includes ACS & tubing as of now */
	  } /* Descent Engine */
	Component ACS {} /* Attitude Control System */
	Component BellyCamera {
	  mass = 2;
	  cost = 0.25;
	}
	Component Radar {
	  mass = 3;
	  cost = 0.5;
	}
	Component Software {}
	Component Disconnect {}
	Component Structure {}
	Component Electronics {}
      } /* propulsion */
  } /* Lander */
  
  Component SRM { /* Solid Rocket Motor */
    transparent mass = sum(key mass of SRM) : 974;
    Cost = 250*old mass*2.2/1E6; /* in millions */
    Component Fuel
      {
	dV = 2560; /* m/s */
	mass = rocket.deltaM(old spacecraft.Mass,
		      290, dV);
      }
    Component drySRM /* Thiokol Star 37FXP */
      {
	mass = 63.6; 
      }
    Component LanderInterface 
      {
	Component Clamp {mass = 0.001*(old icebreaker.mass);}
	Component Harness {mass = 2;}
      }
  } /* SRM */
  } /* Spacecraft */

  Component DSN {cost = 1500/1E6 * duration;}
  
  Component SystemEngineeringDDT {
    overhead = 2.0;
    Component Coreteam {
      personCost = 0.8*Overhead;
      cost = personCost*work; /* 80k / engineer */
      work = 20*3;
    }
    Component facility {cost = 6;}
    Component IronBird {cost = 2;}
    Component prototype {cost = 2;}
    Component mobilityTestbed {cost = 0.6;}
    Component testing {
      Component equipment {cost = 1;}
      Component facility {cost = 0.5;}
      Component radiationFacility {cost = 0.2;}
      Component transportation {cost = 0.6;}
      Component burnInrisk {cost = 3;}
    } /* testing */
    
    Component Management {cost = 0.15*sum(key cost of desc lunarice 
			  except set( management.cost, contingency.cost ) );}
    Component Contingency {cost = 0.2*management.cost/0.15*1.15;}
    
    Component Software {
      unitcost = 0.65*overhead;
      cost = unitcost*sum(key work); /* 60-70 + overhead */
      Component Simulation {work = 5*2;}
      Component CAD {work = 15*2;}
      Component Ground {work = 2*2;}
      Component Science {work = 2*2;}
      Component Test {work = 5*2;}
      Component Injection {work = 5*2;}
      Component Satellite {work = 2*2;}
      Component Lander {work = 5*2;}
      Component HousekeepingSafeguarding {work = 5*2;}
      Component Contingency {work = 5*2;}
    } /* software */
    
  } /* systemEngineeringDDT */
  
} /* lunarice */

/********************** Optional (sub) systems & components ****************/

Component Rocket {
  G = 9.81; /* m/s^2 */
  DeltaM(Mf, Isp, dV) = Mf*(exp(dV/(Isp*earth.G))-1)/exp(dV/(Isp*earth.G));

 /* Returns fuel needed given final Mass after burn, Isp & deltaV
    Do not forget to include rocket dry mass in your Mf, may need 
    iterations, so what does not?

    dV: Delta V 
    Vj: Jet velocity, Vj = Isp*G;
    Mf: Spacecraft mass 
    mu = exp(dV/Vj) 
    Propellant required (approx) = dV/Vj*mo 
    Propellent required (actual) = (mu - 1)/mu * mf 
    Isp: lbf/(lbmass/sec) [in Earth pounds force] 
    */
}

Component OrbiterStage {
  Component Margin {
    mass = 0.2*sum(key mass of OrbiterStage); /* same as power margin */
  }
  
  Component SRM {
    Component Fuel{mass = 33.52;}
    Component rocket {
      mass = 7.5;
      cost = 0.5;
    }
    Component interface {
      mass = 4;
      cost = 0.1;
    }
  } /* SRM */
  Component Satellite {
    Component Solararray {
      mass = 2;
      cost = 0.1;
    }
    Component Gyro {
      mass = 8;
      power = 28;
      cost = 0.5;
    }
    Component computerBoard {
      mass = 2;
      power = 10;
      cost = 0.2;
    }
    Component cameraElectronics {
      mass = 3;
      distrib power = distribution(max=10,min=2);
      cost = 0.4;
    }
    Component CommMoon {
      mass = 3;
      power = 20;
      cost = 0.5;
    }
    Component CommEarth {
      Component horn {
	mass = 1;
      }
      Component comm {
	mass = 5;
	power = 40;
	cost = 1;
      }
      Component gimbal {
	mass = 5;
	power = 5;
	cost = 0.5;
      }
    } /* CommEarth */
    Component battery {
      mass = 5.5;
      cost = 0.2;
    }
    Component PMAD {
      mass = 4;
      cost = 0.3;
    }
    
    Component structure {
      cost = 1;
      Component support {
	mass = 0.2*sum(desc mass of satellite);
      }
      Component bus {
	mass = 5;
      }
      Component tankmount {
	mass = 0.3*sum(desc mass of coldgastank);
      }
    } /* Structure */
    
    Component propulsion {
      cost = 0.5;
      Component coldgastank {
	Component coldgas {mass = 6;}
	Component tank{mass = 2;}
      } /* coldgastank */
      Component nozzles {mass = 2;}
      Component actuators {mass = 2;}
    } /* propulsion */
    
    Component mockup {cost = 0.6;}
  } /* satellite */
} /* orbiter stage */
  


Component CommSatellite{ /* If there is satellite */
  mass = 3;
  power = 10; /* for 1 W output */
  cost = 0.5 + backupEarth.cost;
  Component Transciever{}
  Component Modem{}
  Component Cable{}
  Component Antenna{}
  Component backupEarth{
    mass = 3;
    distrib power = distribution(max=30,avg=0,min=0);
    cost = 1;
  }
}

      





