thomson_problem.cpp

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
using namespace std;

typedef long double ld;
////////////////////////////////////////////////////////////////////////
ld h=0.6;
int n=5;
long long step=100000;
////////////////////////////////////////////////////////////////////////
ld rho_ij(ld theta_i,ld phi_i,ld theta_j,ld phi_j){
    ld tmp=1
        -sin(theta_i)*cos(phi_i)*sin(theta_j)*cos(phi_j)
        -sin(theta_i)*sin(phi_i)*sin(theta_j)*sin(phi_j)
        -cos(theta_i)*cos(theta_j);
    return sqrt(tmp);
}

ld rho_ij_inverse_theta_i(ld theta_i,ld phi_i,ld theta_j,ld phi_j){
    ld tmp=-cos(theta_i)*cos(phi_i)*sin(theta_j)*cos(phi_j)
        -cos(theta_i)*sin(phi_i)*sin(theta_j)*sin(theta_j)
        +sin(theta_i)*cos(theta_j);
    return -tmp/pow(rho_ij(theta_i,phi_i,theta_j,phi_j),3);
}

ld rho_ij_inverse_phi_i(ld theta_i,ld phi_i,ld theta_j,ld phi_j){
    ld tmp=sin(theta_i)*sin(phi_i)*sin(theta_j)*cos(phi_j)
        -sin(theta_i)*cos(phi_i)*sin(theta_j)*sin(phi_j);
    return -tmp/pow(rho_ij(theta_i,phi_i,theta_j,phi_j),3);
}
/////////////////////////////////////////////////////////////////////////
ld potential(vector<vector<ld> > &X){
    ld res=0;
    for(int i=0;i<n;i++){
        for(int j=i+1;j<n;j++){
            res+=1/rho_ij(X[2][i],X[3][i],X[2][j],X[3][j]);
        }
    }
    return res;
}

void ftheta_dd(vector<vector<ld> > &X,vector<ld> &res){
    for(int j=0;j<n;j++){
        res[j]=sin(X[2][j])*cos(X[2][j])*pow(X[1][j],2);
    }
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j) continue;
            for(int k=0;k<n;k++){
                res[k]=rho_ij_inverse_phi_i(X[2][i],X[3][i],X[2][j],X[3][j])/pow(sin(X[2][i]),2);
            }
        }
    }
    return;
}

void fphi_dd(vector<vector<ld> > &X,vector<ld> &res){
    for(int i=0;i<n;i++){
        res[i]=sin(X[2][i])*cos(X[2][i])*pow(X[1][i],2);
        for(int j=0;j<n;j++){
            if(i==j) continue;
            for(int k=0;k<n;k++){
                res[k]-=rho_ij_inverse_phi_i(X[2][i],X[3][i],X[2][j],X[3][j])/pow(sin(X[2][i]),2);
            }
        }
    }
    return;
}
////////////////////////////////////////////////////////////////////////
void function(vector<vector<ld> > &X,vector<vector<ld> > &res){
    ftheta_dd(X,res[0]);
    fphi_dd(X,res[1]);
    for(int i=2;i<=3;i++){
        for(int j=0;j<n;j++){
            res[i][j]=X[i][j];
        }
    }
    return;
}

void add(vector<vector<ld> > &added,vector<vector<ld> > &adder,ld mult){
    for(int i=0;i<4;i++){
        for(int j=0;j<n;j++){
            added[i][j]+=adder[i][j]*mult;
        }
    }
    return;
}
///////////////////////////////////////////////////////////////////////
int main(){
	vector<vector<ld> > prev(4,vector<ld> (n,0.0)),cmp(4,vector<ld> (n,0.0));
    for(int j=0;j<n;j++){
        prev[2][j]=(rand()%RAND_MAX)*M_PI;
    }
    for(int j=0;j<n;j++){
        prev[3][j]=(rand()%RAND_MAX)*2*M_PI;
    }
    
    vector<vector<ld> > ans(4,vector<ld> (n));
    ld potential_min=potential(prev);
    for(long long s=0;s<step;s++){
        vector<vector<ld> > k1(4,vector<ld> (n)),k2(4,vector<ld> (n));
        vector<vector<ld> > k3(4,vector<ld> (n)),k4(4,vector<ld> (n));
        vector<vector<ld> > X1(4,vector<ld> (n,0.0)),X2(4,vector<ld> (n,0.0));
        vector<vector<ld> > X3(4,vector<ld> (n,0.0)),X4(4,vector<ld> (n,0.0));
        X1=X2=X3=X4=prev;
        function(X1,k1);
        add(X2,k1,h/2);
        function(X2,k2);
        add(X3,k2,h/2);
        function(X3,k3);
        add(X4,k3,h);
        function(X4,k4);
        for(int i=0;i<4;i++){
            for(int j=0;j<n;j++){
            cmp[i][j]=prev[i][j]+h/6*(k1[i][j]+2*k2[i][j]+2*k3[i][j]+k4[i][j]);
            }
        }
        for(int j=0;j<n;j++){
            cmp[2][j]=cmp[2][j]-M_PI*floor(cmp[2][j]/M_PI);
        }
        for(int j=0;j<n;j++){
            cmp[3][j]=cmp[3][j]-2*M_PI*floor(cmp[3][j]/2/M_PI);
        }

        ld p=potential(cmp);
        if(p<potential_min){
            potential_min=p;
            ans=cmp;
        }
        /*
        bool ok=true;
        for(int j=0;j<n;j++){
            if(abs(cmp[2][j]-prev[2][j])>pow(h,5) || abs(cmp[3][j]-prev[3][j])>pow(h,5)){
                ok=false;
            }
        }
        if(ok) break;
        */
        if(s%(step/10)==0){
            for(int i=0;i<s/(step/10);i++) cout << "*";
            cout << endl;
            cout << p << endl;
        }
    }
    cout << potential_min << endl;
    for(int j=0;j<n;j++){
        cout << cmp[2][j] << " " << cmp[3][j] << endl;
    }
	return 0;
}