Scusate, ora ho usato il code tag, non so se ora è meglio.
#include <iostream>
#include <math.h>
#include <fstream>
#include <cstdlib>
#include <ctime>
using namespace std;
//parametri
const double mu=398600*pow(10.0,9);//costante gravità terrestre
const double Rt=6378000;//raggio terrestre
const double T=127500;//spinta
const double isp=280;//impulso specifico
const double d=2;//diametro lanciatore
const double rob=7078000;//raggio obiettivo
//condizioni iniziali
const double r0=Rt;//raggio iniziale
const double v0=0;//velocità iniziale
const double gam0=3.1415/2;//angolo di volo iniziale
const double m0=6491;//massa iniziale
const double csi0=0;//anomalia iniziale
const double alpha0=0;//angolo di attacco iniziale
const double g0=9.80665;//accelerazione di gravità sl
//tempo
const double bt=114;//tempo di accensione
const double dt=0.1;//passo di integrazione
//dichiarazione tempo
double tlo,tpo;//variabili da ottimizzare
void traiettoria(double tlo1,double tpo1);//valutazione fitness
double traiettoria1(double tlo,double tpo);//calcolo dati in uscita
void pso();//ottimizzatore
//programma prinicipale
int main(){
pso();
traiettoria(tlo,tpo);
system("pause");
return 0;
}
//ottimizzatore
void pso(){
int ni=200;//numero individui
int np=2;//numero parametri
int nc=1000;//numero cicli
double lbtlo=3;//limite inferiore tempo di lift off
double ubtlo=12;//limite superiore tempo di lift off
double lbtpo=10;//limite inferiore tempo di pitch over
double ubtpo=35;//limite superiore tempo di pitch over
double p[ni][np],v[ni][np],bpg[2],bp[ni][np];
double bfg,fit[ni],bf[ni];
double tlo1,tpo1;
srand(time(0));
for(int i=0;i<ni;i++){
for(int j=0;j<np;j++){
p[i][j]=rand()%1000001/1000000.0;
};
bf[i]=1000;
};
for(int i=0;i<ni;i++){
for(int j=0;j<np;j++){
v[i][j]=rand()%1000001/1000000.0;
};
};
bfg=1000;
for(int j=0;j<nc;j++){
cout <<j<<" "<<bfg<<endl;
for(int i=0;i<ni;i++){
tlo1=(ubtlo-lbtlo)*p[i][1]+lbtlo;
tpo1=(ubtpo-lbtpo)*p[i][2]+lbtpo;
fit[i]=traiettoria1(tlo1,tpo1);
if(fit[i]<bfg){
bfg=fit[i];
bpg[1]=p[i][1];
bpg[2]=p[i][2];
};
if(fit[i]<bf[i]){
bf[i]=fit[i];
bp[i][1]=p[i][1];
bp[i][2]=p[i][2];
};
srand(time(0));
v[i][1]=0.4*v[i][1]+0.5*(rand()%1001 /1000.0)*(bp[i][1]-p[i][1])+0.5*(rand()%1001 /1000.0)*(bpg[1]-p[i][1]);
v[i][1]=0.4*v[i][2]+0.5*(rand()%1001 /1000.0)*(bp[i][2]-p[i][2])+0.5*(rand()%1001 /1000.0)*(bpg[2]-p[i][2]);
p[i][1]=p[i][1]+v[i][1];
p[i][2]=p[i][2]+v[i][2];
};
};
tlo=(ubtlo-lbtlo)*bpg[1]+lbtlo;
tpo=(ubtpo-lbtpo)*bpg[2]+lbtpo;
};
//calcolo fitness
double traiettoria1(double tlo1,double tpo1){
int k,n=bt/dt;
double tspan[n];
double r[n],v[n],gam[n],m[n],csi[n],alpha[n];
double E,h,p,e,a,rmax;
long double fit;
r[0]=r0;
v[0]=v0;
gam[0]=gam0;
m[0]=m0;
csi[0]=csi0;
alpha[0]=alpha0;
tspan[0]=0;
// lift off
for(int i=0;i<tlo1/dt;i++){
r[i+1]=r[i]+v[i]*dt;
v[i+1]=v[i]+(-mu/(pow(r[i],2.0))+T/m[i])*dt;
gam[i+1]=gam[i];
m[i+1]=m[i]-T/(isp*9.80665)*dt;
csi[i+1]=csi[i];
alpha[i+1]=alpha[i];
tspan[i+1]=tspan[i]+dt;
};
// pitch-over
for(int i=tlo1/dt;i<tpo1/dt;i++){
r[i+1]=r[i]+v[i]*sin(gam[i])*dt;
v[i+1]=v[i]+(-mu*sin(gam[i])/(pow(r[i],2.0))+T*cos(alpha[i])/m[i])*dt;
gam[i+1]=gam[i]+((pow(v[i],2.0)/r[i]-mu/pow(r[i],2.0))*cos(gam[i])+T*sin(alpha[i])/m[i])*dt/v[i];
m[i+1]=m[i]-T/(isp*9.80665)*dt;
csi[i+1]=csi[i]+v[i]*cos(gam[i])/r[i]*dt;
alpha[i+1]=alpha[i]+(-3.1415/180+mu*cos(gam[i])/(pow(r[i],2.0)*v[i])-T*sin(alpha[i])/(m[i]*v[i]))*dt;
tspan[i+1]=tspan[i]+dt;
};
//gravity turn
for(int i=tpo1/dt;i<bt/dt;i++){
r[i+1]=r[i]+v[i]*sin(gam[i])*dt;
v[i+1]=v[i]+(-mu*sin(gam[i])/(pow(r[i],2.0))+T/m[i])*dt;
gam[i+1]=gam[i]+((pow(v[i],2.0)/r[i]-mu/pow(r[i],2.0))*cos(gam[i]))*dt/v[i];
m[i+1]=m[i]-T/(isp*9.80665)*dt;
csi[i+1]=csi[i]+v[i]*cos(gam[i])/r[i]*dt;
if(i==tpo/dt){
alpha[i+1]=0;
}
else{
alpha[i+1]=alpha[i];
};
tspan[i+1]=tspan[i]+dt;
};
k=bt/dt;
E=-mu/r[k]+pow(v[k],2.0)/2;
a=-mu/(2*E);
h=r[k]*v[k]*cos(gam[k]);
p=pow(h,2.0)/mu;
e=sqrt(1-p/a);
rmax=a*(1+e);
fit=200*sqrt((rmax-rob,2.0))/rob+E*2*rob/mu;
return fit;
};
void traiettoria(double tlo,double tpo){
int k,n=bt/dt;
double tspan[n];
double r[n],v[n],gam[n],m[n],csi[n],alpha[n];
double E,h,p,e,a,rmax;
ofstream raggio("raggio.txt");
ofstream velo("velocity.txt");
ofstream angolo("angolo di volo.txt");
ofstream massa("massa.txt");
ofstream anomalia("anomalia.txt");
ofstream attacco("angolo di attacco.txt");
r[0]=r0;
v[0]=v0;
gam[0]=gam0;
m[0]=m0;
csi[0]=csi0;
alpha[0]=alpha0;
tspan[0]=0;
raggio<<tspan[0]<<" "<<r[0]-Rt<<endl;
velo<<tspan[0]<<" "<<v[0]<<endl;
angolo<<tspan[0]<<" "<<gam[0]<<endl;
massa<<tspan[0]<<" "<<m[0]<<endl;
anomalia<<tspan[0]<<" "<<csi[0]<<endl;
attacco<<tspan[0]<<" "<<alpha[0]<<endl;
// lift off
for(int i=0;i<tlo/dt;i++){
r[i+1]=r[i]+v[i]*dt;
v[i+1]=v[i]+(-mu/(pow(r[i],2.0))+T/m[i])*dt;
gam[i+1]=gam[i];
m[i+1]=m[i]-T/(isp*9.80665)*dt;
csi[i+1]=csi[i];
alpha[i+1]=alpha[i];
tspan[i+1]=tspan[i]+dt;
raggio<<tspan[i+1]<<" "<<r[i+1]-Rt<<endl;
velo<<tspan[i+1]<<" "<<v[i+1]<<endl;
angolo<<tspan[i+1]<<" "<<gam[i+1]<<endl;
massa<<tspan[i+1]<<" "<<m[i+1]<<endl;
anomalia<<tspan[i+1]<<" "<<csi[i+1]<<endl;
attacco<<tspan[i+1]<<" "<<alpha[i+1]<<endl;
};
// pitch-over
for(int i=tlo/dt;i<tpo/dt;i++){
r[i+1]=r[i]+v[i]*sin(gam[i])*dt;
v[i+1]=v[i]+(-mu*sin(gam[i])/(pow(r[i],2.0))+T*cos(alpha[i])/m[i])*dt;
gam[i+1]=gam[i]+((pow(v[i],2.0)/r[i]-mu/pow(r[i],2.0))*cos(gam[i])+T*sin(alpha[i])/m[i])*dt/v[i];
m[i+1]=m[i]-T/(isp*9.80665)*dt;
csi[i+1]=csi[i]+v[i]*cos(gam[i])/r[i]*dt;
alpha[i+1]=alpha[i]+(-3.1415/180+mu*cos(gam[i])/(pow(r[i],2.0)*v[i])-T*sin(alpha[i])/(m[i]*v[i]))*dt;
tspan[i+1]=tspan[i]+dt;
raggio<<tspan[i+1]<<" "<<r[i+1]-Rt<<endl;
velo<<tspan[i+1]<<" "<<v[i+1]<<endl;
angolo<<tspan[i+1]<<" "<<gam[i+1]<<endl;
massa<<tspan[i+1]<<" "<<m[i+1]<<endl;
anomalia<<tspan[i+1]<<" "<<csi[i+1]<<endl;
attacco<<tspan[i+1]<<" "<<alpha[i+1]<<endl;
};
//gravity turn
for(int i=tpo/dt;i<bt/dt;i++){
r[i+1]=r[i]+v[i]*sin(gam[i])*dt;
v[i+1]=v[i]+(-mu*sin(gam[i])/(pow(r[i],2.0))+T/m[i])*dt;
gam[i+1]=gam[i]+((pow(v[i],2.0)/r[i]-mu/pow(r[i],2.0))*cos(gam[i]))*dt/v[i];
m[i+1]=m[i]-T/(isp*9.80665)*dt;
csi[i+1]=csi[i]+v[i]*cos(gam[i])/r[i]*dt;
if(i==tpo/dt){
alpha[i+1]=0;
}
else{
alpha[i+1]=alpha[i];
};
tspan[i+1]=tspan[i]+dt;
raggio<<tspan[i+1]<<" "<<r[i+1]-Rt<<endl;
velo<<tspan[i+1]<<" "<<v[i+1]<<endl;
angolo<<tspan[i+1]<<" "<<gam[i+1]<<endl;
massa<<tspan[i+1]<<" "<<m[i+1]<<endl;
anomalia<<tspan[i+1]<<" "<<csi[i+1]<<endl;
attacco<<tspan[i+1]<<" "<<alpha[i+1]<<endl;
};
raggio.close();
velo.close();
angolo.close();
massa.close();
anomalia.close();
attacco.close();
k=bt/dt;
E=-mu/r[k]+pow(v[k],2.0)/2;
a=-mu/(2*E);
h=r[k]*v[k]*cos(gam[k]);
p=pow(h,2.0)/mu;
e=sqrt(1-p/a);
rmax=a*(1+e);
cout<<"la quota raggiunta all'apogeo e' "<<rmax-Rt<<endl;
}