Commit 604b2bf0 authored by Daniel Brown's avatar Daniel Brown
Browse files

adding in square aperture analytic knm calculation, set size with r_ap and...

adding in square aperture analytic knm calculation, set size with r_ap and switch to square aperture with conf mirror square_apterture [0/1]
parent 39cc1767
......@@ -418,6 +418,10 @@ typedef struct ABCD {
double D; //!< D element of ABCD matrix
} ABCD_t;
typedef enum aperture_type {
CIRCULAR, SQUARE
} aperture_type_t;
/**
* Structure to define a degree of freedom in a setup
*/
......@@ -734,6 +738,8 @@ typedef struct mirror {
double phi;
double r_aperture; //!< physical size of mirror, infinite if set to 0
aperture_type_t aperture_type;
double x_off;
double y_off;
double Rcx; //!< radius of curvature in x plane (tangential)
......
......@@ -994,6 +994,7 @@ int set_k_mirror(int mirror_index) {
mirror->knm_order[0] = DEFAULT_KNM_ORDER_FIRST;
mirror->knm_order[1] = DEFAULT_KNM_ORDER_SECOND;
mirror->knm_order[2] = DEFAULT_KNM_ORDER_THIRD;
mirror->knm_order[3] = DEFAULT_KNM_ORDER_FOURTH;
}
// the inner product calculations uses the mirrors own knm_q structure
......@@ -1075,6 +1076,13 @@ int set_k_mirror(int mirror_index) {
break;
case 2:
message(" %i Bayer-Helms\n", i+1);
break;
case 3:
if(mirror->aperture_type == CIRCULAR)
message(" %i Circular aperture\n", i+1);
else
message(" %i Square sperture\n", i+1);
break;
case 4:
message(" %i ROMHOM\n", i+1);
......@@ -1101,6 +1109,9 @@ int set_k_mirror(int mirror_index) {
case 2:
message(" - Bayer-Helms\n");
break;
case 3:
message(" - Aperture\n");
break;
case 4:
message(" - ROMHOM\n");
break;
......@@ -1167,10 +1178,14 @@ int set_k_mirror(int mirror_index) {
break;
/* Removing the Aperture analytic K matrix computation
case 3: //aperture
if (!(integrating & 2) && mirror->r_aperture > 0){
fill_mirror_knm_aperture_matrix_analytic(mirror, _nr1);
if(mirror->aperture_type == SQUARE)
fill_mirror_knm_square_aperture(mirror, nr1, nr2);
else
bug_error("not handled");
if (!(mirror->knm_flags & INT_APERTURE) && mirror->map_merged.save_knm_matrices) {
sprintf(buf, "%s_aperture", mirror->map_merged.filename);
......@@ -1178,10 +1193,10 @@ int set_k_mirror(int mirror_index) {
}
// if we have calculated aperture knm analytically do matrix multiplication
mirror_knm_matrix_mult(&(mirror->knm_aperture), &(mirror->knm), &(mirror->knm));
mirror_knm_matrix_mult(&(mirror->knm_aperture), &(mirror->knm_no_rgouy), &(mirror->knm_no_rgouy));
}
break;
*/
case 4:
if (mirror->map_rom) {
......
......@@ -81,6 +81,7 @@
#define DEG 57.2957795130823208767981548142 //!< Convert from rad to deg
#define RAD 0.01745329251994329576923690768485 //!< Convert from deg to rad
#define SQRTTWO 1.41421356237309514547462185874 //!< Definition of sqrt(2)
#define SQRTPI 1.77245385090551602729816748334
#define R2DB 8.685889638065037 // 20 log_10(e^1) Converts the squeezing factor to db
//define sparse matrix solvers
// ddb - values now in solver_t enum
......@@ -372,7 +373,7 @@
// defines number of coupling coefficients matrices there are: maps, bayer-helm and apertures
// Affects both mirrors and beamsplitters
#define NUM_KNM_TYPES 3
#define NUM_KNM_TYPES 4
/* modes for sparese matrix get_..._elements functions */
#define COUNT 0 //!< COUNT mode, prescane elements for CCS storage
......@@ -384,7 +385,8 @@
#define DEFAULT_KNM_ORDER_FIRST 2
#define DEFAULT_KNM_ORDER_SECOND 4
#define DEFAULT_KNM_ORDER_THIRD 1
#define DEFAULT_KNM_ORDER_THIRD 3
#define DEFAULT_KNM_ORDER_FOURTH 1
#endif // KAT_CONSTANTS_H
......
......@@ -13,11 +13,42 @@
#include <stdlib.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
extern init_variables_t init;
extern interferometer_t inter;
static u_nm_accel_t *acc_11_nr1_1, *acc_11_nr1_2;
static u_nm_accel_t *acc_22_nr2_1, *acc_22_nr2_2;
static u_nm_accel_t *acc_21_nr2_1, *acc_21_nr1_2;
static u_nm_accel_t *acc_12_nr1_1, *acc_12_nr2_2;
void alloc_knm_accel_mirror_aperture_mem(long *bytes) {
// need to pick out the maximum number order for allocating memory
int max_m=mem.hg_mode_order, max_n=mem.hg_mode_order;
// make sure they are null before allocating new memory
assert(acc_11_nr1_1 == NULL);
assert(acc_11_nr1_2 == NULL);
assert(acc_12_nr1_1 == NULL);
assert(acc_12_nr2_2 == NULL);
assert(acc_21_nr2_1 == NULL);
assert(acc_21_nr1_2 == NULL);
assert(acc_22_nr2_1 == NULL);
assert(acc_22_nr2_2 == NULL);
// allocate memory for all the accelerators
acc_11_nr1_1 = u_nm_accel_alloc(max_n, max_m, bytes);
acc_11_nr1_2 = u_nm_accel_alloc(max_n, max_m, bytes);
acc_22_nr2_1 = u_nm_accel_alloc(max_n, max_m, bytes);
acc_22_nr2_2 = u_nm_accel_alloc(max_n, max_m, bytes);
acc_12_nr1_1 = u_nm_accel_alloc(max_n, max_m, bytes);
acc_12_nr2_2 = u_nm_accel_alloc(max_n, max_m, bytes);
acc_21_nr2_1 = u_nm_accel_alloc(max_n, max_m, bytes);
acc_21_nr1_2 = u_nm_accel_alloc(max_n, max_m, bytes);
}
double C_n(int n, double x){
assert(n>=0);
......@@ -128,4 +159,107 @@ void fill_mirror_knm_aperture_matrix_analytic(mirror_t *mirror, double nr1){
mirror->knm_aperture.k11[i][j].re = G_matrix_element(m,n,k,l,rhosq11);
}
}
}
\ No newline at end of file
}
complex_t __compute_HG_square(int n, int _n, double R, u_n_accel_t *a1, u_n_accel_t *a2) {
assert(a1 != NULL);
assert(a2 != NULL);
assert(a1->sqrt2_wz == a2->sqrt2_wz);
assert(R > 0);
int p;
complex_t k = z_by_x(z_by_zc(a1->prefac[n], a2->prefac[_n]), 1.0 / a1->sqrt2_wz);
// rescale R for integration limits
R *= a1->sqrt2_wz;
double expR = exp(-R*R);
double sum = 0, tmp = 0;
for (p = 0; p <= min(n, _n); p++) {
int _p = n + _n - 2*p - 1;
if(_p+1 == 0) {
tmp = SQRTPI * erf(R);
} else if((_p+1) % 2 == 1) {
tmp = 0;
} else {
tmp = 2 * hermite(_p, 0) - expR * (hermite(_p, R) - hermite(_p, -R));
}
sum += tmp * pow_two(p) * fac(p) * gsl_sf_choose(n, p) * gsl_sf_choose(_n, p);
}
return z_by_x(k, sum);
}
void fill_mirror_knm_square_aperture(mirror_t *mirror, double nr1, double nr2) {
int max_m, max_n;
get_tem_modes_from_field_index(&max_n, &max_m, inter.num_fields - 1);
max_m = max_n; // max m is no the m is also the max_n. e.g. we get n=1 m=2 for maxtem 2
mirror_knm_q_t *kq = &(mirror->knm_q);
if (CALC_MR_KNM(mirror, 11)) {
u_nm_accel_get(acc_11_nr1_1, max_n, max_m, kq->qxt1_11, kq->qyt1_11, nr1);
u_nm_accel_get(acc_11_nr1_2, max_n, max_m, kq->qxt2_11, kq->qyt2_11, nr1);
}
if (CALC_MR_KNM(mirror, 22)) {
u_nm_accel_get(acc_22_nr2_1, max_n, max_m, kq->qxt1_22, kq->qyt1_22, nr2);
u_nm_accel_get(acc_22_nr2_2, max_n, max_m, kq->qxt2_22, kq->qyt2_22, nr2);
}
if (CALC_MR_KNM(mirror, 12)) {
u_nm_accel_get(acc_12_nr1_1, max_n, max_m, kq->qxt1_12, kq->qyt1_12, nr2);
u_nm_accel_get(acc_12_nr2_2, max_n, max_m, kq->qxt2_12, kq->qyt2_12, nr2);
}
if (CALC_MR_KNM(mirror, 21)) {
u_nm_accel_get(acc_21_nr2_1, max_n, max_m, kq->qxt1_21, kq->qyt1_21, nr1);
u_nm_accel_get(acc_21_nr1_2, max_n, max_m, kq->qxt2_21, kq->qyt2_21, nr1);
}
mirror_knm_t *amap = &(mirror->knm_aperture);
// as we will be calculating aperture knm we need to allocate the memory for it
// remember that free also needs to be called later
if (!IS_MIRROR_KNM_ALLOCD(amap))
bug_error("mirror aperture knm has not been allocated");
mirror->knm_aperture.IsIdentities = false;
int i,j,n,m,k,l;
for(i=0; i<inter.num_fields; i++){
get_tem_modes_from_field_index(&n,&m,i);
for(j=0; j<inter.num_fields; j++){
get_tem_modes_from_field_index(&k,&l,j);
if (CALC_MR_KNM(mirror, 11)) {
mirror->knm_aperture.k11[i][j] = z_by_z(__compute_HG_square(n, k, mirror->r_aperture, acc_11_nr1_1->acc_n, acc_11_nr1_2->acc_n),
__compute_HG_square(m, l, mirror->r_aperture, acc_11_nr1_1->acc_m, acc_11_nr1_2->acc_m));
}
if (CALC_MR_KNM(mirror, 12)) {
mirror->knm_aperture.k12[i][j] = z_by_z(__compute_HG_square(n, k, mirror->r_aperture, acc_12_nr1_1->acc_n, acc_12_nr2_2->acc_n),
__compute_HG_square(m, l, mirror->r_aperture, acc_12_nr1_1->acc_m, acc_12_nr2_2->acc_m));
}
if (CALC_MR_KNM(mirror, 21)) {
mirror->knm_aperture.k21[i][j] = z_by_z(__compute_HG_square(n, k, mirror->r_aperture, acc_21_nr2_1->acc_n, acc_21_nr1_2->acc_n),
__compute_HG_square(m, l, mirror->r_aperture, acc_21_nr2_1->acc_m, acc_21_nr1_2->acc_m));
}
if (CALC_MR_KNM(mirror, 22)) {
mirror->knm_aperture.k22[i][j] = z_by_z(__compute_HG_square(n, k, mirror->r_aperture, acc_22_nr2_1->acc_n, acc_22_nr2_2->acc_n),
__compute_HG_square(m, l, mirror->r_aperture, acc_22_nr2_1->acc_m, acc_22_nr2_2->acc_m));
}
}
}
}
......@@ -8,7 +8,10 @@
#ifndef KAT_KNM_APERTURE_H
#define KAT_KNM_APERTURE_H
void alloc_knm_accel_mirror_aperture_mem(long *bytes);
void fill_mirror_knm_aperture_matrix_analytic(mirror_t *mirror, double nr1);
void fill_mirror_knm_square_aperture(mirror_t *mirror, double nr1, double nr2);
#endif /* KAT_KNM_APERTURE_H */
......@@ -320,8 +320,8 @@ bool should_integrate_mirror_aperture_knm(mirror_t *mirror, int mismatch, int as
} else {
// otherwise we need to calculate the aperture coupling coefficient matrix here analytically
if(!(mirror->map_merged.usermessages & APERTURE_NO_ANALYTIC)){
warn("Analytic aperture coupling coefficient calculation not completely implemented yet, becareful!\n");
if(!(mirror->map_merged.usermessages & APERTURE_NO_ANALYTIC )&& mirror->aperture_type == CIRCULAR){
warn("Analytic circular aperture coupling coefficient calculation not completely implemented yet, becareful!\n");
mirror->map_merged.usermessages |= APERTURE_NO_ANALYTIC;
}
......
......@@ -204,7 +204,8 @@ void alloc_matrix_ccs_nnz(long *bytes, ifo_matrix_vars_t *matrix_ccs){
gerror("Could not allocate CCS row_idx");
} else
*bytes += try_bytes;
#if INCLUDE_PARALUTION == 1
if(options.sparse_solver == PARALUTION) {
// Paralution needs COO format to assemble matrix
if ((matrix_ccs->M.col_idx = calloc(NNZ, size)) == NULL) {
......@@ -214,6 +215,7 @@ void alloc_matrix_ccs_nnz(long *bytes, ifo_matrix_vars_t *matrix_ccs){
} else {
matrix_ccs->M.col_idx = NULL;
}
#endif
try_bytes = 2 * NNZ * sizeof (double);
......@@ -4257,7 +4259,8 @@ void build_ccs_matrix(ifo_matrix_vars_t *matrix){
// by definition the last eqn + 1 must be the total number of nnz
set_col_ptr(matrix->M.col_ptr, eqns, matrix->M.num_nonzeros);
#if INCLUDE_PARALUTION == 1
// Now build the column index storage needed by certain solvers
if (options.sparse_solver == PARALUTION) {
int column = 0, i = 0, tmp = 0;
......@@ -4271,6 +4274,7 @@ void build_ccs_matrix(ifo_matrix_vars_t *matrix){
tmp += matrix->M.row_col_count[column];
}
}
#endif
}
void get_ccs_matrix(ifo_matrix_vars_t *matrix){
......@@ -4585,7 +4589,7 @@ void solve_ccs_matrix(ifo_matrix_vars_t *matrix, double *rhs, bool transpose, bo
handle_NICSLU(NicsLUc_SolveFast(&objs->nicslu, (complex__t*)(void*)rhs));
#endif
#if INCLUDE_NICSLU == 1
#if INCLUDE_PARALUTION == 1
} else if(options.sparse_solver == PARALUTION) {
solve_paralution_matrix(matrix, rhs, transpose, conjugate);
#endif
......@@ -4608,10 +4612,14 @@ void solve_ccs_matrix_sparse_input(ifo_matrix_vars_t *matrix, double *rhs, bool
if(ko->Common.status != 0)
gerror("KLU error whilst solving: %i", ko->Common.status);
#if INCLUDE_NICSLU == 1
} else if(options.sparse_solver == NICSLU){
gerror("NICSLU doesn't support conjugate solving at the moment...");
#endif
#if INCLUDE_PARALUTION == 1
} else if(options.sparse_solver == PARALUTION){
bug_error("NOT DONE");
#endif
} else
bug_error("Sparse solver not handled");
}
......
......@@ -37,6 +37,7 @@
#include "kat_dump.h"
#include "kat_aa.h"
#include "kat_knm_mirror.h"
#include "kat_knm_aperture.h"
#include "kat_knm_bs.h"
#include "kat_quant.h"
......@@ -1409,6 +1410,7 @@ int allocate_memory_for_mirror_list(long int *bytes) {
*bytes += (mem.num_mirrors + 1) * sizeof (mirror_t);
alloc_knm_accel_mirror_aperture_mem(bytes);
alloc_knm_accel_mirror_mem(bytes);
mirror_knm_alloc(&mrtmap, bytes);
......
......@@ -5,6 +5,7 @@
* Created on 10 March 2015, 16:06
*/
#if INCLUDE_PARALUTION == 1
#ifndef KAT_PARALUTION_H
#define KAT_PARALUTION_H
......@@ -25,4 +26,5 @@ extern "C" {
}
#endif
#endif /* KAT_PARALUTION_H */
#endif
......@@ -1597,6 +1597,7 @@ void read_mirror(const char *command_string, int mode) {
mirror->knm_order[0] = DEFAULT_KNM_ORDER_FIRST;
mirror->knm_order[1] = DEFAULT_KNM_ORDER_SECOND;
mirror->knm_order[2] = DEFAULT_KNM_ORDER_THIRD;
mirror->knm_order[3] = DEFAULT_KNM_ORDER_FOURTH;
// default to changing n and q over the bayer-helms computation
mirror->knm_change_q = DEFAULT_KNM_CHANGE_Q;
......@@ -1696,6 +1697,7 @@ void read_mirror2(const char *command_string) {
mirror->knm_order[0] = DEFAULT_KNM_ORDER_FIRST;
mirror->knm_order[1] = DEFAULT_KNM_ORDER_SECOND;
mirror->knm_order[2] = DEFAULT_KNM_ORDER_THIRD;
mirror->knm_order[3] = DEFAULT_KNM_ORDER_FOURTH;
mirror->x_off = 0.0;
mirror->y_off = 0.0;
......@@ -1844,6 +1846,7 @@ void read_beamsplitter(const char *command_string, int mode) {
bs->knm_order[0] = DEFAULT_KNM_ORDER_FIRST;
bs->knm_order[1] = DEFAULT_KNM_ORDER_SECOND;
bs->knm_order[2] = DEFAULT_KNM_ORDER_THIRD;
bs->knm_order[3] = DEFAULT_KNM_ORDER_FOURTH;
bs->knm_change_q = DEFAULT_KNM_CHANGE_Q;
// default to integrating aperture until analytic version put in
......@@ -1967,6 +1970,7 @@ void read_beamsplitter2(const char *command_string) {
bs->knm_order[0] = DEFAULT_KNM_ORDER_FIRST;
bs->knm_order[1] = DEFAULT_KNM_ORDER_SECOND;
bs->knm_order[2] = DEFAULT_KNM_ORDER_THIRD;
bs->knm_order[3] = DEFAULT_KNM_ORDER_FOURTH;
bs->knm_change_q = DEFAULT_KNM_CHANGE_Q;
// default to integrating aperture until analytic version put in
......@@ -6971,6 +6975,17 @@ void read_conf(const char *command_string) {
gerror("knm_force_saved must be either 1 or 0 for true or false\n");
mirror->knm_force_saved = (bool) ival;
} else if (strcmp(config_name, "square_aperture") == 0) {
GET_INT
if (ival != 1 && ival != 0)
gerror("square_aperture must be either 1 or 0 for true or false\n");
if (ival == 0)
mirror->aperture_type = CIRCULAR;
else
mirror->aperture_type = SQUARE;
} else if (strcmp(config_name, "knm_cuba_cores") == 0) {
GET_INT
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment